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Axisymmetric,  thin-layer  Navier-Stokes  equations  with 
radiative  heat  transfer  and  internal  heat  generation  terms 
are  used  to  describe  the  thermal  and  flow  fields  of  interest 
to  Ultrahigh  Temperature  Vapor  Core  Reactor  systems.  An 
implicit-explicit,  finite  volume  method  in  conjunction  with 
an  algebraic  two-layer  eddy  viscosity  turbulence  model  is 
used  to  numerically  solve  these  equations.  Algebraic  grid 
clustering  technique  is  utilized  to  resolve  the  thin 
boundary  layer  near  the  solid  wall.  Convective  heat  fluxes 
into  the  solid  wall  are  predicted  by  Fourier's  law  in  the 
highly  stretched  grid  systems.  Radiative  heat  transport 
process  is  modeled  using  both  the  1— D integral  approximation 
and  the  Rosseland' s diffusion  approximation.  Using  a 
straight  tube  geometry  and  the  converging— diverging  geometry 
of  the  reactor  exit  nozzle,  the  convective  and  radiative 
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heat  transfer  rates  are  calculated.  Existing  experimental 
data  for  high  temperature  air  flow  in  tubes  are  used  to 
benchmark  the  model.  Good  agreement  between  the  numerical 
results  and  the  measured  data  indicates  the  validity  of  the 
developed  method  for  the  intended  analysis. 

Results  of  the  analysis  indicate  that  the  accurate 
prediction  of  convective  heat  transfer  rates  at  the  wall 
requires  grid  sizes  less  than  the  thickness  of  the  laminar 
sublayer.  For  the  converging-diverging  nozzle,  the  maximum 
value  of  convective  and  radiative  heat  fluxes  occurs  just 
upstream  to  the  throat  of  the  nozzle  where  the  mass  flux  is 
maximum.  In  general,  the  convective  heat  flux  increases 
with  increasing  the  stagnation  pressure  which  corresponds  to 
larger  mass  fluxes.  It  is  also  found  that  the  diffusion 
approximation  for  the  calculation  of  radiative  heat  flux  is 
only  valid  for  an  optically  thick  gas.  For  less  opaque  gas 
conditions  the  integral  transport  approximation  is  applied. 
The  radiative  heat  transfer  rate  is  dominant  at  low 
pressures  and  high  gas  temperatures. 
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CHAPTER  1 
INTRODUCTION 

1 . 1 Motivation  and  Background 

Of  the  numerous  space  nuclear  power  concepts  proposed 
and  studied  in  recent  years,  a UF4-fuel  based  vapor  core 
reactor  (UTVR)  with  magnetohydrodynamic  (MHD)  generator  is 
one  of  the  most  promising  and  interesting.  The  reactor 
concept  is  employed  for  the  analysis  of  fluid  flow  and  heat 
transfer  problems  in  this  study.  The  analysis  is  especially 
focused  on  parts  of  the  UTVR  system  where  the  maximum  heat 
loading  to  the  wall  is  expected. 

Figure  1.1  shows  the  basic  concept  of  ultrahigh 
temperature  vapor  reactor  with  advanced  energy  conversion 
system.  In  this  system,  the  average  exit  temperature  range 
for  the  core  is  between  4000-5000  K and  in  the  MHD  duct  is 
between  2100-2500  K.  The  inner  wall  temperatures  of  the 
core  and  the  MHD  generator  are  around  2000  K [1]. 

Since  a very  high  temperature  working  fluid  and  fuel 
are  used  to  enhance  the  efficiency  of  the  energy  conversion 
system,  potential  structural  and  thermal  design  problems 
arise  due  to  the  high  heat  transfer  expected  to  occur 
between  the  gas  and  the  solid  walls  of  the  core  and  the 
nozzle.  Accordingly,  the  accurate  prediction  of  local  wall 
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Figure  1.1  200  MWe  Ultrahigh  Temperature  Vapor  Reactor-MHD  Generator  Rankine  Space 

Power  System 
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heat  transfer  rates  is  of  crucial  importance  in  order  to 
prevent  a catastrophic  failure,  such  as  a deterioration  of 
wall  structural  integrity.  The  energy  deposition  to  the 
solid  wall  can  be  made  up  of  fission  fragment  heating,  p and 
Y ray  heating,  radiative  heat  transfer,  and  convective  wall 
heat  transfer.  However,  only  the  convective  and  radiative 
heat  transfer  mechanisms,  which  dominate  the  other  heat 
transfer  mechanisms,  are  considered  in  this  study. 

Prediction  of  the  convective  heat  transfer  rate  in 
complex  flow  fields  similar  to  those  of  the  UTVR-MHD  system 
is  cumbersome.  Moreover,  the  radiative  heat  transfer 
process  due  to  a high  temperature  flowing  medium  is  another 
energy  transfer  mechanism  which  adds  to  the  complexity  of 
the  problem.  Consequently,  a consistent  method  of 
predicting  the  convective  and  radiative  heat  transfer  rates 
must  be  developed.  This  gives  a strong  motivation  for  the 
present  work. 

Prediction  of  convective  heat  transfer  rates  has  been  a 
major  area  of  research  for  several  decades  in  heat  transfer. 
Initially,  dimensional  analysis  was  used  to  develop 
empirical  correlations  for  heat  transfer  coefficients. 
However,  these  empirical  correlations  are  limited  by  the 
system  geometry  and  the  fully  developed  flow  condition. 

For  a fully-developed  turbulent  flow  in  a tube,  Dittus- 
Boelter  [2],  Colburn  [3],  and  Sieder  and  Tate  [4]  developed 
a variety  of  empirical  correlations.  These  investigations 
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primarily  used  the  data  from  experiments  performed  in  long 
tubes  under  low  temperature  flow  conditions  at  the  ambient 
pressure.  Later,  Zellnik  and  Churchill  [5]  extended  the 
temperature  range  of  the  flowing  medium,  and  partially 
considered  the  entrance  effect.  However,  the  limitations  of 
the  above  methods  still  exist  for  geometries  and  flow 
conditions . 

With  the  advent  of  fast  computers,  computational  fluid 
flow  and  heat  transfer  models,  which  can  be  applied  to  every 
possible  flow  geometry  and  condition,  have  been  developed 
and  refined.  However,  there  are  still  many  problems  to  be 
solved  for  the  prediction  of  wall  heat  transfer  by 
computational  heat  transfer  modeling. 

In  particular,  previous  computational  fluid  flow  and 
heat  transfer  modeling  for  the  wall  heat  transfer  prediction 
have  focused  on  the  question  of  providing  detailed  near-wall 
resolution.  Viegas,  Rubesin,  and  Horstman  [6]  have 
addressed  the  difficulty  of  dealing  with  the  numerical 
computation  of  wall  heat  transfer  prediction.  Thus,  they 
have  employed  the  empirical  or  semiempirical  information 
concerning  the  nature  of  the  turbulent  boundary  layer  near  a 
surface  instead  using  a fine  near— wall  grid. 

Other  researchers  used  the  empirical  or  semi-empirical 
information,  that  is,  "wall  functions"  for  their  heat 
transfer  studies  [7-10].  Rose  [11]  also  used  the  wall 
function  in  conjunction  with  an  algebraic  turbulence  model 
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in  his  heat  transfer  modeling.  However,  a recent  report 
[12]  pointed  out  the  disadvantages  of  the  wall  function 
method,  which  are  too  sensitive  to  the  near  wall  grid 
position,  and  the  inaccuracy  at  the  region  of  flow 
separation  and  reattachment. 

In  this  report  a transfer  model  is  presented  based  upon 
a numerical  patching  technique  in  which  two  different  sets 
of  governing  equations  for  a boundary  layer  region  and  an 
outer  flow  region  are  set  up,  respectively.  In  the  patching 
method,  the  boundary  layer  solution  is  coupled  to  the  outer 
flow  solution  through  a matching  condition  applied  at  the 
outer  edge  of  the  wall  layer.  This  method  improves 
computational  efficiency  by  using  a coarse  grid  in  the  outer 
layer  region  and  a very  fine  grid  in  the  wall  layer  region. 

As  a result,  wall  heat  flux  and  wall  shear  stress  are 
calculated  directly  from  the  slopes  of  the  calculated 
velocity  and  temperature  profiles  at  the  wall.  However,  the 
patching  technique  is  not  appropriate  for  the  nozzle  flow 
which  is  characterized  by  strong  viscous/inviscid, 
shock/boundary  layer  interactions. 

It  is,  therefore,  required  to  use  an  approach  which 
solves  a suitable  subset  of  Navier-Stokes  equations  since 
this  method  considers  these  interactions  in  a fully  coupled 
manner.  To  use  this  approach  for  prediction  of  heat 
transfer  rate  to  the  wall,  a highly— stretched  fine  grid  near 
the  solid  wall  is  required.  The  fine  mesh,  however,  yields 
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computational  difficulties  in  terms  of  the  stability 
limitation  and  the  computational  time.  It  is  logical  then 
to  explore  numerical  methods  which  would  relieve  the 
severities  of  these  limitations. 

The  numerical  solution  which  requires  fine  spatial  mesh 
distributions  to  resolve  viscous  layer  near-wall  boundaries 
is  subject  to  the  small  time  step  limitation  for  the 
explicit  algorithm.  To  remove  the  time  step  limitation,  the 
fully  implicit  numerical  schemes  were  developed  in  the  mid- 
1970's  by  Briley  and  McDonald  [13]  and  Beam-Warming  [14]. 
Although  the  implicit  approximate  factorization  procedure 
can  have  a larger  time  step  than  the  explicit  procedure,  the 
factorization  introduces  an  error  proportional  to  (At)2 
which  restricts  the  time  step  from  being  too  large. 

An  hybrid  implicit-explicit  method  presented  by 
McCormack  [15,16]  requires  the  simple  inversion  of  block 
bidiagonal  matrices  and,  hence,  greatly  reduces  the 
computing  time  needed.  This  numerical  scheme  with  Baldwin 
and  Lomax's  algebraic  turbulence  model  [17]  is  utilized  in 
this  study. 

In  addition,  radiative  heat  transfer  caused  by  the  high 
temperature  fuel-working  fluid  becomes  another  important 
transport  mechanism  in  this  system.  In  the  feasibility 
study  of  the  direct  flow  gaseous  core  reactor  system  [18],  a 
diffusion  approximation  for  an  optically  dense  gas  appeared 
to  be  a good  approximation  as  well  as  an  integration  method 
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of  the  radiation  equation.  Other  researchers  [19-22]  also 
used  the  diffusion  approximation  and  mentioned  that  the 
method  was  valid  in  the  optically  thick  medium,  that  is,  the 
mean  penetration  distance  of  photon  is  quite  small  compared 
to  the  characteristic  dimension  of  the  medium. 

With  the  above  considerations,  the  model,  which  solves 
a suitable  subset  of  Navier-Stokes  equations  with  radiative 
heat  transfer  and  an  internal  heat  source  term  in  a very 
fine  near-wall  grid,  is  regarded  as  one  of  the  best  models 
for  this  system.  This  information,  based  on  the  previous 
studies,  gives  the  direction  of  modeling  to  fit  the  present 
flow  phenomena  and  problem  characteristics.  In  the  following 
section,  the  objectives  of  the  present  work  and  a preview  of 
the  mathematical  model  will  be  described  in  more  detail. 

1 . 2 Objectives  and  Preview  of  Mathematical  Model 

The  primary  objective  of  this  work  is  to  develop  a 
computational  fluid  dynamics  and  heat  transfer  model  for 
flow  of  a very  high  temperature  radiating  gas  with  internal 
heat  generation  in  a UTVR/MHD  system.  In  particular,  to 
provide  a methodology  for  accurate  estimation  of  heat 
transfer  rates  from  the  hot  gas  to  the  inner  wall  of  the 
reactor  chamber  and  associated  nozzle  is  a primary  concern 
of  this  study. 

To  achieve  this  goal,  a two-dimensional  (2-D) , 
axisymmetric  Navier-Stokes  equations  solver  along  with  the 
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turbulence  model  is  used.  A near-wall  diffusion  model  for 
convective  wall  heat  transfer  and  Rosseland' s diffusion 
model  and  1-D  integral  approximation  for  radiative  heat 
transfer  are  also  employed. 

The  procedure  is  based  on  the  solution  of  Navier-Stokes 
equations  for  unsteady,  compressible,  viscous,  turbulent 
flow.  To  obtain  the  flow  and  thermal  field  solution,  the 
axisymmetric  thin  layer  Navier-Stokes  equations  with 
radiative  heat  transfer  and  internal  heat  source  are  used 
and  numerically  solved. 

A hybrid  implicit-explicit  method  in  conjunction  with 
Lomax  and  Baldwin' s two-layer  algebraic  turbulence  model  is 
used,  and  the  finite  volume  approach  is  also  applied  to  this 
analysis.  A flux-splitting  procedure  [14],  which  takes 
advantage  of  the  natural  direction  of  information  travel  and 
enables  the  implicit  coefficient  matrix  to  be  more 
diagonally  dominant,  is  employed.  The  Rosseland' s diffusion 
approximation  is  employed  to  account  for  the  radiation 
contribution,  which  models  the  radiation  heat  transfer  in 
the  partial  differential  equation,  instead  of  the  partial 
integro-dif ferential  equation.  The  1-D  integral 
approximation  is  also  used  for  the  radiative  heat  transfer. 

For  prediction  of  convective  heat  transfer  rate,  a 
highly  stretched  grid  system,  in  which  several  grid  points 
can  be  located  in  the  viscous  sublayer,  is  used.  Thus,  the 
convective  heat  transfer  rate  to  the  wall  is  directly 
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determined  by  the  known  temperature  gradient  at  the  wall. 
With  the  developed  model,  the  effect  of  internal  heat 
generation  on  the  fluid  flow  and  heat  transfer  in  a straight 
tube  is  investigated.  The  effect  of  numerical  grid 
characteristics,  i.e.,  skewness  and  orthogonality,  is  also 
examined  in  the  flow  and  thermal  field  solutions.  In 
addition,  more  computations  are  carried  out  for  a variety  of 
flow  conditions,  geometric  parameters,  radiative  mean 
opacity,  and  constant  power  densities. 

1 . 3 Overview 

Theoretical  modeling,  including  the  thin  layer  Navier- 
Stokes  equation,  diffusion  model  for  radiative  heat  transfer 
and  turbulence  model,  is  described  in  Chapter  2.  Chapter  3 
provides  an  algebraic  stretching  grid  scheme  for  a 
converging-diverging  nozzle  and  also  explains  the  algebraic 
method  for  obtaining  an  orthogonal  grid  system  in  the 
nozzle.  Then  in  Chapter  4 a numerical  algorithm  and 
boundary  conditions  are  described  for  the  solution  of  the 
governing  equations.  Chapter  5 provides  a method  for 
evaluating  the  heat  transfer  using  the  known  temperature 
profile.  Finally,  Chapter  6 presents  results  and  Chapter  7 
provides  conclusions. 


CHAPTER  2 

MATHEMATICAL  FORMULATION 

In  this  chapter  the  governing  equations  for  unsteady 
compressible  flow  with  thermal  radiation  and  internal  heat 
generation  are  summarized.  Based  on  the  physical 
considerations  as  well  as  computer  limitations,  the  thin- 
layer  Navier-Stokes  equations  are  used  for  this  study.  The 
two-dimensional  (2-D)  axisymmetric  thin-layer  Navier-Stokes 
equations  are  expressed  in  conservation  law  form.  The 
diffusion  model  for  radiant  heat  transfer  and  the  Baldwin- 
Lomax  turbulence  model  are  described. 

2 . 1 Governing  Equation 

In  modeling  the  simultaneous  development  of  momentum 
and  thermal  boundary  layer  of  a radiating  gas  in  the 
axisymmetric  cylindrical  coordinate  system  (z,r),  the 
following  assumptions  are  made: 

(1)  Viscous  dissipation,  heat  conduction,  and  thermal 
radiation  in  the  streamwise  direction  are 
negligible . 

(2)  Optical  properties  of  the  participating  gas  are 
approximated  in  terms  of  those  of  a grey  emitting, 
absorbing,  and  nonscattering  gas. 
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(3)  There  is  no  line  radiation  at  temperatures 
considered  for  this  analysis. 

(4)  When  the  gas  properties  at  high  temperature  are 
not  available,  the  perfect  gas  law  is  applied  to 
evaluate  them. 

(5)  There  is  no  dissociation  or  phase  change  for  UF4 
gas  in  the  1800  to  4000  K temperature  range.  This 
assumption  is  based  on  the  following  experimental 
data  taken  from  Reference  [23] . 

(a)  Dissociation  of  UF„  is  less  than  1%  at  1 atm 
and  3500  K. 

(b)  Dissociation  of  UF4  is  less  than  1%  at  100  atm 
and  5000  K. 

Under  these  assumptions,  the  differential  equations  used  to 
describe  the  mean  flow  for  this  study  are  the  time- 
dependent,  mass-averaged  Navier-Stokes  equations  with 
radiative  and  constant  internal  heat  source  terms  for  the 
axially  symmetric  flow  of  a compressible  fluid.  Utilizing 
the  Boussinesq  assumption  for  the  turbulence  model,  the 
resulting  equations  [24,  25],  in  strong  conservation  form, 
can  be  written  as  follows: 


du 

dt 


♦ Ei 

dz 


(2.1) 


in  which 
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and  the  viscous  terms  Gv  and  Hv  are 
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Here,  p is  the  density;  u and  v are  the  velocity  components 
in  the  z and  r directions,  respectively;  Q is  the  uniform 
internal  heat  generation  term;  q"c  and  q"r 


are  conductive 
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and  radiative  heat  fluxes,  respectively;  p is  the  pressure; 
and  e is  the  total  energy  per  unit  volume  which  is  related 
to  the  internal  energy  £ and  the  kinetic  energy.  The  total 
energy  per  unit  volume  is  written  by 


The  total  viscosity  |iT  is  given  as  the  sum  of  its 
molecular,  ji2,  and  eddy  components,  pt: 


The  molecular  viscosity  for  air  is  assumed  to  obey 
Sutherland's  formula  [25]: 


where  Cj  and  C2  are  constants  for  a given  gas.  For  air  at 
moderate  temperature,  q = 1 . 458x10**  {kg/m  s y/K)  and 

C2  = 110.4  K.  The  turbulent  eddy  viscosity  is  obtained  from 

the  turbulent  eddy  viscosity  model  proposed  by  Baldwin  and 
Lomax.  The  heat  flux  terms  are  comprised  of  the  conductive 

heat  flux  term,  q" , and  the  radiative  heat  flux  term,  q'J . 


e = p [ e + (u2  + v2  ) ] 


(2.5) 


Hr  = Hi  + Ht 


(2.6) 


(2.7) 


The  conductive  heat  flux  can  be  written  as 
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dy 


(2.8) 


where 


(2.9) 


Typical  values  of  the  molecular  and  turbulent  Prandtl 
numbers  for  air  at  standard  conditions  are  0.9  and  0.72, 
respectively.  The  radiative  heat  flux  term  is  presented  in 
the  following  section. 

The  complete  set  of  time-dependent  thin-layer  Navier- 
Stokes  equations  (2.1)  includes  one  continuity  equation,  two 
momentum  equations,  and  one  energy  equation,  but  there  are 
six  unknowns,  p,  u,  v,  e,  p,  and  T.  In  order  to  close  the 
system  of  equations,  the  fluid  is  assumed  to  be  a perfect 
gas.  The  equation  of  state  and  other  relations  for  a 
perfect  gas  are 


where  R is  the  gas  constant,  h is  enthalpy,  and  cp  and  cv 
are  the  specific  heat  at  constant  pressure  and  volume, 
respectively.  y is  the  ratio  of  specific  heats,  which  for 
air  at  standard  conditions  is  1.4.  Using  the  above  formula, 
(2.5)  and  (2.10),  two  additional  equations  relating  p and  T 


p = pRT,  e = cvT,  h = cpT,  y = 


(2.10) 
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are  obtained  as  follows: 


P = ( y —1 ) [e  - p (u2  + v2)  ] 

(2.11) 

T = (?  ' 1}  [-  ~ \ (u2  + v2)] 

R p 2 

Hence,  equations  (2.1)  and  (2.11)  are  the  governing 
equations  for  axisymmetric,  unsteady  compressible  flow 
problems.  In  this  study,  two  sets  of  fluid  properties  are 
used.  One  is  air  at  standard  conditions  for  the 
verification  purpose  of  fluid  flow  and  heat  transfer 
solutions,  and  the  other  one  is  UF4  gas  properties  evaluated 
at  2000  K.  The  properties  of  the  UF4  gas  are  presented  in 
the  following  subsection. 

2.1.1  UF„  Properties 

As  mentioned  above,  two  sets  of  fluid  properties  are 
utilized  to  analyze  the  fluid  flow  and  heat  transfer 
problems  in  this  work.  One  is  the  air  and  the  other  is  pure 
UF4  gas  which  is  described  in  this  section.  Pure  UF4  gas  is 
treated  as  a perfect  gas  and  is  assumed  not  to  be 
dissociated  at  temperatures  between  1800  K and  4000  K.  The 
properties  of  pure  UF4  gas  are  summarized  in  Table  2.1  and 
Table  2.2.  The  specific  heats  and  the  ratio  of  specific 
heats  (y)  data  are  taken  from  Reference  [26]  . The  Viscosity 
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2. 

1 Pure  UF4  properties  at  2000  K. 

Mw 

(Kg/Kg-  mol) 

311.15 

R 

(J/Kg-K) 

26.72 

y 

1.08 

cP 

(J/Kg-K) 

345.95 

Cv 

(J/Kg-K) 

319.31 

Hi 

(Kg/m-  s) 

8.67  x 

Ki 

(W/m-  K) 

0.028 

Pr 

0.915 

Table  2.2  Heat  capacities  of  UF4  gas  at  temperature  between 
1800  and  4000  K. 


Temperature 
in  K 

Cp  in  J 
Kg'1-  K'1 

Cv  in  J 
Kg'1-  K'1 

Ratio  of 

Specific  Heat,  y 

1800 

345.68 

318.90 

1.08 

2000 

345.95 

319.31 

1.08 

2200 

346.22 

319.58 

1.08 

2400 

346.29 

319.71 

1.08 

2600 

346.62 

319.98 

1.08 

2800 

346.76 

320.12 

1.08 

3000 

346.89 

320.16 

1.08 

3200 

347.03 

320.25 

1.08 

3400 

347.03 

320.25 

1.08 

3600 

347.16 

320.38 

1.08 

3800 

347.16 

320.38 

1.08 

4000 

347.16 

320.52 

1.08 
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of  pure  UF4  gas  is  assumed  to  be  the  same  as  that  of  UF6 
reported  by  Oliver  and  Dugan  [27] . All  properties  are 
estimated  at  2000  K.  Indeed,  it  matters  little  in  this  work 
what  particular  gas  properties  are  employed,  because  the 
primary  concern  of  this  study  is  not  in  the  application  of  a 
design  problem  but  in  providing  a robust  and  accurate  heat 
transfer  model.  Using  the  known  ratio  of  specific  heats,  y, 
a laminar  Prandtl  number  is  determined  by  Euken's  formula 
[28] 


Pr 


4Y 

9y  - 5 


(2.12) 


The  Prandtl  number  is  often  used  to  determine  the 
coefficient  of  thermal  conductivity  kc  once  the  viscosity  is 
known. 

2 . 2 Radiative  Heat  Transfer 

The  thin-layer  Navier-Stokes  equations  in  conjunction 
with  the  Baldwin-Lomax  algebraic  turbulence  model  were 
expressed  and  used  to  solve  flow  equations  in  a coupled 
manner  with  the  energy  equation.  Internal  heat  generation 
and  radiative  heat  transfer  are  considered  in  the  energy 

equation  as  source  terms,  Q and  V-g",  respectively.  In  the 

present  study,  the  internal  heat  source  term  is  assumed  to 
be  constant;  however,  modeling  or  some  simplification  for 
radiative  heat  transfer  is  needed.  The  assumption  of  opaque 
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gas  can  lead  to  a substantial  simplification  such  as  a 
diffusion  approximation.  In  this  diffusion  approximation, 
the  radiation  arriving  at  any  location  comes  only  from  the 
immediate  surroundings  because  any  other  radiation  is 
absorbed  before  arriving  at  that  location.  The  diffusion 
approximation  with  the  above  assumption  leads  to  a 
considerable  simplification  in  the  expression  for  radiative 
heat  flux,  but  the  simplicity  is  offset  by  its  disadvantage. 
For  instance,  the  approximation  is  valid  at  points  far  from 
the  solid  surface  and  only  for  intense  absorption,  i.e., 
when  the  fractional  variation  in  temperature  is  small  in  a 
distance  of  one  mean  free  path.  Additionally,  the 
approximation  breaks  down  in  the  vicinity  of  the  surface, 
since  it  does  not  take  into  account  radiation  leaving  from 
the  surface.  Therefore,  the  limit  and  the  accuracy  of  the 
approximation  should  be  examined.  For  comparison  purposes, 
the  radiative  heat  flux  is  predicted  by  two  different 
approximations,  the  diffusion  approximation  and  the  one- 
dimensional ( 1 — D ) equation  of  radiative  transfer.  The  two 
approximations  are  described  below. 

2.2.1  Diffusion  Approximation 

As  mentioned  above,  the  diffusion  approximation  is  only 
valid  for  an  optically  dense  medium.  This  can  be  stated 
more  rigorously  by  letting  H be  a path  length  over  which  the 
radiant  energy  density  does  change  appreciably,  and  letting 
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lm  be  the  extinction  mean  free  path.  Then,  the  diffusion 
approximation  is  applicable  for  the  condition  of  lm/H  « 1. 
Under  the  assumption  of  an  optically  thick  medium,  the 
radiative  heat  flux  is  proportional  to  the  temperature 
gradient  and  is  written  as  [29] 


__16oT^Vr  _ _ktjt 
3 aR 


(2.13) 


where  Kr,  the  radiative  conductivity,  is  defined  by 


Kw 


16  or* 

3aR 


(2.14) 


where  aR  is  the  Rosseland  mean  opacity  and  a is  the  Stefan- 

Boltzmann  constant.  To  obtain  the  mean  opacity,  the 
spectral  absorption  coefficient  must  be  averaged  over  the 
spectrum.  Consequently,  the  Rosseland  mean  opacity  of  a 
nonscattering  gas  is  defined  by  [29] 


a* 


L 


- dBv 
~dr 


cfv 


/« 


- n 2 


o a„  dT 


dv 


(2.15) 


where  v is  the  photon  wave  number,  Bv  is  the  Planck 
function,  av  is  the  spectral  absorption  coefficient,  and  n 


20 


is  the  real  index  of  refraction.  Due  to  the  lack  of 
experimental  data  for  the  spectral  absorption  coefficient  of 
UF4,  direct  calculation  of  the  Rosseland  mean  opacity  for 
UF4  is  not  feasible.  However,  for  this  study  where  the 
self-absorption  by  UF4  is  under  consideration,  the  Rosseland 
mean  opacity  can  be  estimated  by  [30] 

aR  ~ NOph  (2.16) 

where  N is  the  molecular  number  density  of  the  gas  and  aph 
is  the  photon  collision  cross  section  per  molecule.  Using 
the  assumption  of  the  perfect  gas  law,  the  molecular  number 
density  can  be  computed  by 


N = 


P_ 
K T 


(2.17) 


where  k is  Boltzmann's  constant,  p is  the  gas  pressure,  and 
T is  the  gas  temperature.  Therefore,  the  radiative 
conductivity  can  be  written  as 


kr 


16osbk 

3<JPtP 


4 


(2.18) 


The  photon  absorption  cross  section  for  UF4  gas  is  very 
close  to  that  of  UF6  and  approximated  based  on  the  data 
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presented  in  Reference  31.  The  frequency  averaged  photon 
absorption  cross  section  in  the  2000  A to  4200  A wavelength 
range  for  UF4  is  estimated  to  be  equal  to  the  frequency 
averaged  photon  absorption  cross  section  of  UF6  which  is 
2.76  million  barns  (2.76  x 10"22  m2)  . The  frequency  averaged 
Rosseland  opacity  for  UF4  gas  based  on  the  estimated  photon 
absorption  cross  section  in  the  temperature  range  of  2000  K 
to  5000  K is 


aR  = 2x10 4 - 

A (J7 


(2.19) 


where  p is  in  atm,  T is  in  K,  and  aR  is  in  cm-1.  The 
estimated  values  of  the  frequency  averaged  Rosseland  opacity 
agree  with  the  published  values  of  the  mean  opacity  for  pure 
uranium  gas  at  5000  K [32].  Using  the  estimated  opacities, 
the  Rosseland  approximation  leads  to  a considerable 
simplification  in  the  expression  for  radiative  heat  flux. 
However,  when  the  frequency  averaged  Rosseland  coefficient 
is  applied  to  real  gases  which  are  usually  transparent  in 
some  wavelength  regions,  some  errors  might  be  introduced. 

To  assess  the  sensitivity  of  radiative  transfer  to  the 
average  value  of  the  Rosseland  absorption  coefficient,  a 
wavelength-band  application  is  performed  for  two  wavelength- 
bands,  2000  A to  3000  A and  3000  A to  4200  A.  in  addition, 
the  diffusion  approximation  is  not  accurate  near  the  solid 


boundary  since  it  does  not  take  into  account  radiation 
leaving  from  the  surface;  the  presence  of  the  surface 
changes  the  coefficient  from  16/3  to  8/3  [33], 
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2.2.2  Approximation  by  Using  1-D  Equation  of  Radiative 
Transfer 

This  approximation  is  performed  with  the  integral  form 
of  the  one-dimensional  equation  of  radiative  transfer  by 
using  the  frequency  averaged  opacity  utilized  in  the 
diffusion  approximation.  With  no  incident  radiation  at  the 
center  line,  the  equation  of  radiative  transfer  is  written 
as  [29] 


where  i is  the  intensity,  and  a is  the  frequency  averaged 
Rosseland  opacity,  and  I is  the  source  function,  which  can 
be  reduced  to  the  local  blackbody  intensity  in  absorbing 
medium  without  scattering.  The  local  blackbody  intensity  is 
given  by 


In  order  to  simplify  the  equation  (2.20),  moreover,  it  is 
assumed  that  the  wall  is  a blackbody,  and  the  radiation 
crossing  a unit  area  normal  to  the  r direction  is  only 


(2.20) 


or4 


(2.21) 
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considered  in  intensity  terms  of  the  positive  and  negative 
sense.  The  net  radiative  flux  in  the  gas  in  the  positive  r 
direction  is 

g/  = <7+  - gr-  = 71  (i  + (r)  - i-(r))  (2.22) 

where  q+  is  the  radiative  flux  traveling  from  gas  to  the 
boundary  and  q~  is  the  radiative  flux  in  the  opposite 
direction . 

The  assumption  of  considering  only  the  radiation  normal  to 
the  z direction  might  possibly  be  acceptable  under  certain 
conditions,  that  is,  the  temperature  variation  in  the  axial 
direction  is  very  small  and  the  geometry  is  very  simple, 
such  as  a straight  tube. 

The  simplified  integral  equation  shown  above,  is  solved 
by  using  the  simple  trapezoidal  rule.  At  every  node  point 
the  net  radiative  flux  is  calculated  by  summing  the  positive 
values  evaluated  from  the  interlayer  of  the  gas  to  the  wall 
and  the  negative  values  for  the  opposite  direction.  In  this 
summation,  the  radiation  is  considered  to  be  absorbed  within 
a five-times  distance  of  the  mean  free  path. 

2 . 3 Algebraic  Turbulence  Model 

The  algebraic  turbulence  model  employed  in  this  study 
is  the  one  proposed  by  Baldwin  and  Lomax  [17]  . It  is  a two- 
layer  algebraic  eddy  viscosity  model  in  which  is 
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calculated  for  an  inner  and  an  outer  region.  The  inner 
region  follows  the  Prandtl-Van  Driest  formulation.  In  both 
the  inner  and  outer  formulation,  the  distribution  of 
vorticity  is  used  to  determine  the  length  scales,  thereby 
avoiding  the  necessity  of  finding  the  outer  edge  of  the 
boundary  layer.  For  the  inner  region, 

- P-12M  (2-23) 


where 


1 = Ky  [1  - exp  ( -y*/A*)  ] 


(2.24) 


and 


y*  = p^y  = /w 


(2.25) 


Here,  y is  the  normal  distance  to  the  surface  and  |u|  is  the 
absolute  magnitude  of  vorticity.  The  nondimensional 
parameter,  y* , is  the  law  of  the  wall  coordinate.  1 is  the 
mixing  length  and  uT,  the  frictional  velocity,  is 


u 


T 


(2.26) 
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where  xw,  pw  and  are  the  local  shear  stress,  density,  and 

laminar  viscosity  evaluated  at  the  wall. 

The  eddy  viscosity  for  the  outer  region  is  given  by 


( outer  KCcpPFwakoFklel)(y')  (2.27) 

where  K is  the  Clauser  constant,  Ccp  is  an  additional 
constant,  and  Fwalce  is  the  smaller  value  of  the  following: 


wake 


F wake  y max^max 

(~'wkyaaxudif  /-^max  • 


(2.28) 


The  quantities  yraax  and  Fmax  are  determined  from  the  function 


-F(y)  =y|u|[l  - exp  {-y*/A*)  ] . (2.29) 

In  a wake,  the  exponential  term  of  the  above  equation  is  set 
equal  to  zero.  The  quantity  Fmax  is  the  maximum  value  of 
F (y)  that  occurs  in  a profile,  and  ymax  is  the  value  of  y at 
which  it  occurs.  The  function  Fkleb{y)  is  the  Klebanoff 
intermittency  factor  given  by 

Fklot  (y>  = [1  + 5.5  (_Cl-iab7)6]  -1  . 

y max 


(2.30) 
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The  quantity  Udif  is  given  as  the  difference  between  maximum 
and  minimum  total  velocity  in  the  profile 


Udif  = - (v/u^v2)^  . 


(2.31) 


The  minimum  total  velocity  is  taken  to  be  zero,  unless  a 
wake  is  considered.  The  constants  appearing  in  the 
foregoing  relations  have  been  determined  by  requiring 
agreement  with  the  Cebeci  model  for  constant  pressure 
boundary  layers  at  transonic  speeds.  The  constants  used  for 
this  model  are  A+  = 26,  Ccp  =1.6,  Cwk  = 0.25,  Ckleb  = 0.3,  K = 
0.4,  and  K = 0.0168. 


CHAPTER  3 
GRID  GENERATION 

3 . 1 Grid  for  Numerical  Solution 

To  solve  the  governing  equations  in  differential  forms 
for  fluid  mechanics  and  heat  transfer,  an  approximation  to 
the  partial  differential  equations  is  applied  and  a 
numerical  solution  scheme  is  introduced.  The  approximation 
converts  the  partial  derivatives  to  finite  difference 
expressions,  which  are  used  to  rewrite  the  partial 
differential  equations  as  algebraic  equations.  Then  the 
approximate  algebraic  equations,  referred  to  as  finite 
difference  equations,  are  subsequently  solved  at  discrete 
points  within  the  domain  of  interest.  Therefore,  a set  of 
grid  points  within  the  domain,  as  well  as  the  boundaries  of 
the  domain,  must  be  specified.  The  creation  of  such  a grid 
system  is  known  as  grid  generation. 

The  generation  of  a grid  is  one  of  the  central  problems 
in  computing  numerical  solutions.  A well  constructed  grid 
greatly  simplifies  the  numerical  solution  of  systems  of 
partial  differential  equations.  On  the  other  hand,  an 
improper  grid  choice  may  lead  to  instabilities  and  lack  of 
convergence . 
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In  general,  grid  generation  techniques  may  be 
classified  as 

(1)  algebraic  methods, 

(2)  partial  differential  methods, 

(3)  complex  variable  methods  (conformal  mapping) . 

The  selection  of  the  most  suitable  and  efficient  technique 
of  grid  generation  for  a particular  problem  depends  on  the 
geometry  of  the  physical  domain. 

In  the  axisymmetric  code  used  for  this  study,  the 
algebraic  grid  generation  technique  has  been  adopted.  In 
general,  algebraic  methods  use  some  special  formula  for  the 
smooth  grid  network  of  a certain  geometry.  The  converging- 
diverging  nozzle,  which  is  a critical  component  in  the  UTVR- 
MHD  power  system  concept,  with  large  momentum  and  energy 
transfer  gradients  is  considered  in  this  study.  The 
algebraic  grid  generation  technique  for  this  geometry  is 
briefly  explained  in  the  following. 

The  general  shape  of  a converging-diverging  nozzle  is 
shown  in  Figure  3.1.  First,  the  coordinate  of  every  grid 
point  between  the  center  line  and  the  solid  wall  at  the 
throat  is  calculated  by  the  following  formula: 


r{j)  =-i?(  l- 


exp  [ ckb  (j-2)  ] 
Bk 


j-2, ,Jmax  (3.1) 


where  the  parameters,  Ck  and  Ekj  which  are  determined  in 
iterative  procedures,  control  the  degree  of  the  density  or 
the  sparsity  of  the  grid  system.  Another  parameter,  8,  is 
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the  uniform  grid  interval  at  the  throat  for  which  the  throat 
radius  R is  divided  by  the  number  of  radial  direction  grids. 
Then,  the  coordinates  of  grid  points  on  the  solid  wall  are 
algebraically  determined  by  the  uphill  angle,  0,  the 
downhill  angle,  e,  and  the  curvature,  R' , at  the  throat. 

The  coordinate  of  each  grid  point  in  the  radial  direction  is 
determined  by  the  distance  ratio  of  every  grid  interval 
between  the  consecutive  grids  to  the  throat  radius  at  the 
throat.  In  the  same  manner,  all  the  interior  grid 
coordinates  along  the  flow  stream  direction  are  calculated. 

on  the  other  hand,  the  physical  coordinate  grid  must  be 
finer  near  the  walls  and  through  the  regions  of  major 
disturbances,  such  as  the  expansion  region.  Well  away  from 
these  regions  a coarser  physical  grid  can  be  used  for 
computational  efficiency,  while  still  adequately  resolving 
the  flow  field.  The  grid  points  along  the  radial  direction 
are  clustered  toward  the  rigid  flow  boundaries,  where  large 
gradients  are  expected,  using  exponential  stretching.  The 
degree  of  stretching  along  the  radial  direction  are 
controlled  by  the  parameters,  Ck and  Ek,  used  in  Eq.  (3.1). 

The  parameters  are  determined  by  the  analytical  formulas 


Fk-dy-R 


exp  (Ck6)  -l 
exp  ( Ck)  -l 


(3.2) 


and 
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(exp  (Cfc)  -l)  6 exp  ( ck)  - (exp  (Ck6)  -1 ) exp  ( Ck) 


(3.3) 


(exp  (C*.)  -l)2 


Ck  is  the  last  updated  value  for  the  specified  number  of 
iteration  as  follows: 


Ek  is  evaluated  using  the  determined  Ck  in  the  following: 


Attention  must  be  paid  to  the  manner  in  which  the  grid 
points  are  distributed  along  the  axis  because  the 
distribution  also  affects  to  a large  extent  the  accuracy  of 
the  calculations.  The  grid  points  along  the  axis  must  be 
concentrated  in  the  areas  of  large  wall  slopes  and  axial 
locations  where  separation  is  expected  to  occur.  The  same 
grid-clustering  scheme  is  applied  to  the  axial  direction  as 
used  in  the  radial  direction.  The  physical  domain  grid 
constructed  using  Eqs . (3.1)  , (3.2),  (3.3),  (3.4),  and 

(3.5)  is  shown  in  Figure  3.2. 


(3.4) 


£^=exp  ( Ck)  -1.0 


(3.5) 
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Figure  3.2  A skewed  grid  in  a converging-diverging  nozzle 


Figure  3.3  A skewed  grid  in  a bell-type  nozzle. 


33 


The  above  basic  algebraic  technique  can  also  be  applied 
to  a different  type  of  nozzle,  such  as  a bell-type  nozzle, 
with  minor  modification.  The  uniform  grid  system  for  the 
bell-type  nozzle  is  presented  in  Figure  3.3.  Both  types  of 
nozzles  are  used  in  this  work  for  the  study  of  flow  and  heat 
transfer  in  the  nozzle. 

3 . 2 Orthogonal  Grid  with  an  Algebraic  Method 

The  algebraic  method  for  grid  generation  is  getting 
popular  due  to  its  simplicity  and  the  computational  speed  it 
provides.  However,  the  algebraic  grid  system  usually  has 
significant  skewness  which  is  the  main  disadvantage  of  the 
algebraic  method.  Severe  skewness  must  be  avoided  because 
it  creates  problems  with  the  numerical  implementation  when  a 
generalized  coordinates  transformation  is  used  to  transform 
the  physical  domain  grid  to  the  computational  domain  grid. 

In  the  following  a simple  algebraic  method  is  presented 
to  transform  a skewed  grid  to  an  orthogonal  grid.  The 
method  is  applied  to  the  converging-diverging  nozzle.  The 
algebraic  grid  generation  technique  explained  in  the 
previous  section  does  a piecewise  linear  approximation  of 
the  boundary  curve,  S^,  as  shown  in  Figure  3.4.  From  each 
axial  position  along  the  boundary  curve,  a vertical  line  to 
the  symmetry  axis  is  drawn.  Along  each  vertical  line  the 
9-^id  points  are  distributed  closer  to  the  boundary  curve 
using  the  exponential  clustering  scheme. 
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Figure  3.4  Skewed  grid(dotted  line)  and  corrected  grid(solid  line) 
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Figure  3.4  shows  that  this  approach  results  in  severely 
skewed  grids  when  the  slope  of  the  boundary  curve  with 
respect  to  the  symmetry  axis  is  large.  In  order  to  correct 
the  grid  skewness,  it  is  necessary  to  draw  a curve  which  is 
normal  on  both  the  boundary  curve  and  the  axis  of  symmetry. 

For  the  particular  case  where  the  symmetry  axis  is  a 
straight  line  the  curve  is  a section  of  a circle.  The 
center  lies  at  the  intersection  of  the  axis  with  the 
straight  line  segment,  B^,  approximating  the  curved 
boundary.  The  radius  of  the  circle  is  0Blf  which  is  the 
distance  between  the  interaction  point  0 of  B,B2  and  the 
symmetry  axis.  This  circular  sector,  indeed,  intersects  at 
right  angles  with  both  the  boundary  curve  and  the  axis. 

The  grid  points  on  the  vertical  line  presented  in 
Figure  3.5  are  now  shifted  to  new  positions  on  the  circular 
arc.  The  new  positions  are  the  intersections  of  the 
circular  arc  and  the  family  of  lines  emanating  from  the 
center  0 and  passing  through  the  old  positions  on  the 
vertical.  The  radius  of  the  circle  and  the  new  shifted 
coordinates  of  the  grid  points  are  calculated  from  the 
following  trigonometric  relations: 


cc=tan-1  [ {BxB[-B2Bi)  /&&] 
OB^B^Bl—1 


smo 
P =tan'1  (P^l/OBl) 


(3.6) 
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Figure  3.5  Displacement  of  grid  points  from  the  old  position  PI  to 
the  new  position  P'l  (corrected  grid) 
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The  modified  orthogonal  grids  for  the  above  two  types 
of  nozzle  are  shown  in  Figures  3.6  and  3.7,  respectively. 

The  improvement  in  terms  of  orthogonality  is  observed  from  a 
visual  comparison  of  the  skewed  and  the  orthogonal  grid  or 
from  a comparison  of  the  numerical  values  of  the 
interaction  angles  of  the  skewed  and  normalized  grids  given 
in  Table  3.1. 
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Figure  3.6 


Corrected  grid  in  a converging-diverging  nozzle. 


Figure  3.7  Corrected  grid  in  a bell-type  nozzle. 
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Table  3.1  Intersection  angles  of  the  skewed  grid  and  the 

normalized  grid  in  the  nozzle  of  Figures  3.2  and  3.6, 
along  a mesh  line  close  the  curved  rigid  boundary 
where  the  maximum  skewness  is  observed. 


Skewed  Grid 

Orthogonal  Grid 

60° 

90° 

60° 

90° 

60° 

90° 

60° 

90° 

65° 

90° 

73° 

90° 

78° 

90° 

83° 

90° 

85° 

90° 

88° 

90° 

90° 

90° 

89° 

90° 

88° 

90° 

85° 

90° 

83° 

90° 

80° 

90° 

77° 

90° 

76° 

90° 

75° 

90° 

75° 

90° 

75° 

90° 

75° 

90° 

CHAPTER  4 

NUMERICAL  ALGORITHMS 

In  this  chapter  the  finite  volume  discretization 
procedure  which  has  been  followed  to  obtain  a suitable 
algorithm  for  computer  implementation,  is  presented.  A 
hybrid  implicit-explicit  numerical  method  based  on  a flux 
splitting  scheme  is  described.  Finally,  boundary  and 
initial  conditions  and  the  convergence  criterion  are 
discussed. 


4 . 1 Finite  Volume  Discretization 

Finite  volume  methods  approximate  integral  forms  of 
conservation  equations  on  finite  cells.  The  equations  are 
approximated  by  summing  the  fluxes  of  mass,  momentum,  and 
energy  from  neighboring  cells  into  each  cell,  thereby 
satisfying  the  conservation  equations  over  that  cell.  In 
this  method,  the  integral  formulation  of  the  conservation 
laws  are  discretized  directly  in  the  physical  domain  which 
allows  the  use  of  the  shock-capturing  technique  in  the 
transonic  and  supersonic  flow  regions. 

Another  feature  of  the  finite  volume  approach  is  the 
decoupling  of  the  discretization  of  the  equations  from  the 
grid  generation  so  that  the  rigorous,  smooth  variation  of 
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mesh  size  and  an  explicit  coordinate  transformation  are  not 
required  [34].  In  order  to  help  the  understanding  of  the 
process,  consider  an  element  located  at  (i,j)  with  its  four 
boundaries  ab,  be,  cd,  and  da,  as  shown  in  Figure  4.1.  The 
integral  form  of  the  governing  equation  given  in  equation 
(2.1)  can  be  written  as 


and  where  0 is  the  azimuthal  angle  and  A represents  the  area 
which  is  positive  for  counter-clockwise  elements  and 
negative  for  clockwise  elements.  Here,  consider  the  left- 
hand  side,  i.e.,  the  inviscid  terms,  such  that  the  equation 
is  reduced  to  the  Euler  equation.  Application  of  Gauss's 
theorem  yields 


(4.1) 


where 


dv±  j = {r±  jA0)  dA1  j 


+ = 0 (4.2) 


To  express  the  above  integral  form  to  a discrete  form,  let 
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• 
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• 

i,  2 

• 

i+1,2 

/ 
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• 

i-1- 1 

• 

i,  1 

• 

i+1 , 1 

Figure  4.1  A finite  volume  mesh. 
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us  consider  a simple  control  volume  described  in  Figure  4.1 
and  a cell  (i,j).  Then,  equation  (4.2)  can  be  reduced  to 
the  discrete  form: 


En  Ar  + D±  Fa  Az)  + H±  j At  (4.3) 

A±.l 

where  D±[E)  and  D±(F)  i ^ are  the  increments  of  inviscid 
flux  terms  which  can  be  approximated  by  forward  difference 


DAE)itj  = (E)i+1<J  - 

D.<ft,.,  - (ft,.,.!  - (ft,.. 


and  backward  difference 


DAE)  ± j = ( E)itj  - (E)  i_lfj 


(4.4) 


(4.5) 


The  superscript  n to  indicate  present  time  level  is  dropped 
for  convenience. 

4 . 2 Implicit-Explicit  Scheme 


The  hybrid  implicit-explicit  method  is  employed  to 
solve  the  governing  equation  given  in  equation  (2.1).  The 
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basic  algorithm  of  the  method  is  based  on  the  explicit 
predictor-corrector  scheme  proposed  by  MacCormack  [35] . The 
explicit  formulation  is 


A Ufa  = -At 


Vi  + 5(Gi  - Gv) 


A z 


A r 


- H 


(4.6) 


i.J 


where  H = Hv  - Hi  and  Dt  represents  either  forward  or 

backward  and  the  arrows  indicating  vector  quantities  are 
omitted.  The  above  explicit  algorithm  is  numerically  stable 
for  Courant-Friedrich-Lewy  (CFL)  numbers  less  than  or  equal 
to  one.  To  remove  the  CFL  restriction,  an  implicit 
procedure  is  required. 

Modern  implicit  numerical  methods  are  of  the  form  [16] 


[NUMERICS]  6 U±'j  = [PHYSICS] 

The  above  formula  calculates  the  change  in  the  solution 
® ul*j  at  grid  point  i,j  during  the  interval  from  time  nAt  to 

time  (n+l)At.  The  term  on  the  right-hand  side  (physics) 
represents  a finite  difference  or  volume  approximation  to 
the  complete  governing  equation  (2.1)  using  the  solution  at 
time  nAt,  i.e., 
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[PHYSICS]  = ACT/j 
= -A  t 


D±Fd  + D±  ( G±  - Gv)  _ ^ 


A z 


A r 


To  derive  an  implicit  form,  equation  (2.1)  is  differentiated 
with  respect  to  time  as 


»<£> 

3t 


3z 


3r 


dH 

3t 


(4.7) 


where 


dF, 

■W'  s< 


dG^ 

"0C7 


and 


5. 


3Gv 

0(7 


are  the  Jacobians  of 


Fif  Gx  and  Gv  with  respect  to  U. 
Letting 


At 


At  (|£)n  = 
3t 


A Hn 


and 


6(7n+1, 


an  implicit  difference  approximation  to  equation  (4.7)  is 
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J + At 


D • Ai 
A z 


+ At 


D • (Bi  - Bv) 


A r 


Ji.J 


tuft 


= hUf'j  + 


where  the  dots  indicate  that  the  difference  operator  D 
applies  to  all  factors  to  the  right  of  the  operator,  i.e., 

to  the  factor  aoEf. 

To  improve  numerical  efficiency  and  to  take  advantage 
of  natural  directions  in  which  information  propagates,  the 
Steger  and  Warming's  flux  vector  splitting  technique  [36]  is 
used.  Since  the  inviscid  flux  terms  which  represent  the 
macromotion  of  working  fluid  are  homogeneous  of  degree  one 
in  the  element  of  U,  Steger  and  Warming  pointed  out  that 


F1  = AjU 
G±  = B±U 


The  Jacobian  matrices  AL,  can  be  diagonalized  by 
similarity  transformations  Sx,  Sy  such  that 
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A,  = Sx1AaSx 
Bi  = 


where 


i 

0 

0 

-1/c2' 

' 1 

0 

0 

0‘ 

J = 

0 

pc  0 

1 

-u/p 

1/P 

0 

0 

X 

0 

0 

1 

0 

-v/p 

0 

1/p 

0 

p 

-pc  0 

1 

. 

-up 

-vP 

p 

1 

0 

0 

-l/c2 

' 1 

0 

0 

O' 

= 

0 

1 

0 

0 

-u/p 

1/p 

0 

0 

y 

0 

0 pc 

1 

-v/p 

0 

1/p 

0 

0 

0 -pc 

1 

ca. 

s 

-up 

-vp 

p 

and 


u 

0 

0 

0 

V 

0 

0 

0 ‘ 

0 

u+c  0 

0 

0 

V 

0 

0 

0 

0 

u 

0 

tP" 

ii 

0 

0 

v+c 

0 

0 

0 

0 

u-c 

0 

0 

0 

v-c 

Here,  c is  the  speed  of  sound,  a=  (u2  + v2)/2,  p=y-i,  and  y 
is  the  ratio  of  the  specific  heats  of  gas.  The  A,  and  A„ 

A B 

are  diagonal  matrices,  and  the  diagonal  elements  in  these 
matrices  are  the  eigenvalues  of  the  original  Jacobian 
Matrices  B^.  In  general,  some  of  the  elements  of  the 

diagonalized  matrix  are  positive  and  others  are  negative. 
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Their  signs  determine  the  direction  of  information  travel. 
According  to  the  sign  of  these  eigenvalues,  Ai  and  can  be 
written  as 


Ai  ~ Ai+  + 

Bi  = Blt.  + 


where  A1+  and  Bi+  include  only  positive  diagonal  elements  and 
Ai-  and  Bj_  include  the  negative  ones.  Thus,  the  original 
matrices  are  split  as 


Ai  + Bx  A A+Sx 
Ai-  = Sx  A A_SX 
Bi  + ~ <Sy  A 3+<Sy 
5i-  = 5;1AB_sy 


As  a result,  the  fluxes  passing  through  a cell  surface,  for 
instance,  at  i+l/2,j,  can  be  evaluated  as  follows  (see 
Fig. 4 . 2) , 


l+ 1/2.  j Ai+Ui'j  + Ai_Ui+lij 

1+ 1/2,  j = Bl*Ul.j  + Bl-Ul*i.j 


Considering  the  flux  vector  splitting  technique,  the  hybrid 
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Figure  4.2  Direction  of  information  travel  in  time 
to  surface  i+l/2,j. 
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Figure  4.3  Sweep  direction  for  line  Gauss- 
Seidel  iteration. 
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implicit-explicit  algorithm  is  finally  summarized  as 

Predictor : 


[I  + A t( 


= -At 

DSi 

Az 

>'A-  + 

£_  * A+ 

Az 

Az 

At 

fir 

(SGV>  ] 

Ar 


- H 


i.J 


) + At( 


D.  * B.  D.  • B. 


Ar 


Ar 


) (4.8) 


Ctf.?  = Oh  + Sctf.* 


Corrector: 


Attfj  = -At 


D_  • Fv 


Az 


+ S.  • (Gi  - Gy) 


Ar 


[j  + At(^V^i  + + At(£l7l5:  + £^l£l)  (4.9) 


Az 
At 


Az 

n* i * rTn+l 


Ar  Ar 

-J7  »t?.?  - Aof?  * aSi", 

t£?  - ♦ oS7  * »o?;,1) 


where 

(1)  the  superscripts  n,  n+T,  and  n+1  refer  to 
present,  predicted,  and  new  solution  values, 

(2)  the  subscripts  i and  j refer  to  a special  set  of 
mesh  point  locations,  (iAx,  jAy) , 

(3)  Au  is  the  temporal  change  in  the  solution  during 


time  interval  At  and  is  evaluated  by  a local 
explicit  finite  difference  approximation, 
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(4)  the  difference  operators,  £>_•  and  X , 

represent  forward,  backward,  and  central 
difference  approximations  applied  to  all  factors 
to  the  right  of  the  operator  with  respect  to  the 
direction  indicated  in  the  denominator  of  the 
difference  operator  quotient,  i.e., 


The  implicit  equation  of  each  step  can  be  written  as 


where  the  coefficients  A,  B,  C,  D,  and  E are  4x4  block 
matrices.  Equation  (4.10)  can  be  solved  by  the  line  Gauss- 
Seidel  iterative  procedure  [37,38].  The  stencil  of  mesh 
points  used  is  shown  in  Figure  4.3.  A line  at  a time  is 
solved  for  simultaneously  by  using  the  latest  available 
values  ahead  of  and  behind  the  line.  The  parameter  k is  the 
iteration  index.  The  relaxation  algorithm  is  implemented  by 


+ is  equivalent  to 


(4.10) 
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alternately  sweeping  in  the  backward  and  forward  streamwise 
directions.  Alternating  the  direction  of  the  implicit  line 
inversion  has  been  proven  to  be  beneficial  to  the  solution 
of  the  Navier-Stokes  equations  with  highly  stretched  grids. 
It  is  also  consistent  with  the  nature  of  transonic  and 
viscous  separated  flow  [39] . 

4 . 3 Boundary  and  Initial  Conditions 

For  finite  volume  schemes,  the  reflection  condition  has 
often  been  used  to  treat  the  flow  on  a wall  boundary.  The 
boundary  conditions  are  summarized  as 


U, 


i,l 


'i.  1 


Li,l 


Ui.2 

slip 

~Ui.2 

no  slip 

(2^21  ~ 

^.2 

moving  wall 

[~Vi.2 

no  slip 

- 

Vi.2 

moving  wall] 

{ Ti.2 

adiabatic  1 

l2  Tva  11 

Ti.2 

iso  thermal j 

(4.11) 


where  U and  V are  the  velocities  of  the  z and  r directions, 
respectively.  The  constant  stagnation  pressures  and 
temperatures  at  the  inflow  and  outflow  boundary  are  used. 
Initially,  the  flow  is  at  rest  with  pressure  and  temperature 
set  everywhere  to  their  total  values  at  the  entrance.  The 
following  medium  starts  to  move  by  setting  the  exit  pressure 
to  a lower  pressure  than  the  inlet  pressure. 
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4 . 4 Convergence  Criterion 

Since  an  iterative  solution  procedure  is  used  to  solve 
the  equations,  it  is  necessary  to  establish  a convergence 
criterion  which  measures  the  degree  to  which  a computed 
solution  satisfies  the  finite  difference  equation.  In  the 
present  work,  this  convergence  criterion  is  based  on  the 
value  of  velocity  residual  defined  by 


imax  jmax 

E E <[<,  - u?;,1] 2 * t - v#)2)1'2 

R = .'•i  

( imax ) ( jmax ) 


(4.12) 


When  the  reduction  of  the  velocity  residual  is  saturated  to 
an  asymptotic  value,  the  solution  is  considered  converged. 
In  general,  the  asymptotic  value  of  the  residual  varies 
according  to  the  flow  conditions  and  the  geometries  used, 
and  the  values  obtained  in  this  study  were  mostly  less  than 
the  value  of  R = 3 x 10"3.  The  convergence  was  also  able  to 
be  graphically  determined. 


CHAPTER  5 

EVALUATION  OF  HEAT  TRANSFER  RATES 
5 . 1 Convective  Heat  Transfer  Rate 

There  are  two  ways  to  determine  the  convective  heat 
flux  into  the  wall  boundary.  One  way  is  to  use  Fourier's 
law  as  the  product  of  the  fluid  conductivity  and  temperature 
gradient  near  the  wall.  The  other  way  is  to  use  the 
convective  heat  transfer  coefficients  which  can  be  obtained 
from  the  literature  for  fully  developed  flows  of  various 
types  in  conduits  or  in  simple  geometries.  For  the  first 
method,  a very  fine  grid  near  the  wall  is  required  so  that 
the  viscous  sublayer  can  be  properly  resolved  and  then  the 
heat  flux  is  directly  determined  as 


Q"  = <5.1> 

where  the  thermal  conductivity,  kc,  is  evaluated  at  the 
average  temperature  of  the  cell  to  the  wall.  The 
temperature  gradient  in  equation  (5.1)  is  numerically 
calculated  using  the  estimated  cell  temperature  and  the  wall 
temperature.  Usually  the  convective  heat  transfer  is 
predicted  with  the  convective  heat  transfer  coefficient,  hc, 
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as  defined  by 


h 


C 


(5.2) 


where  Tw  is  the  wall  temperature  and  Tb  is  the  bulk 
temperature.  The  bulk  temperature  is  averaged  over  the  flow 
cross-sectional  area  according  to  the  following  relation 


where  A is  the  cross-sectional  flow  area,  Cp  is  the  specific 
heat  at  constant  pressure,  p is  the  density,  and  u is  the 
axial  flow  velocity.  On  the  other  hand,  the  empirical 
correlations  can  be  used  to  calculate  the  convective  heat 
transfer  coefficients.  A typical  correlation  for  the  wall 
convective  heat  transfer  coefficient  for  fully  developed 
turbulent  flow  in  pipes  is  due  to  Dittus  and  Boelter  [2] 


j jPCpUTdA 


(5.3) 


N = 0.023  i?e°-8Pr 1/3 


(5.4) 


where  the  Nusselt  number  is  defined  as  hD/kc.  The  Nusselt 
number  obtained  with  the  above  equation  (5.4)  is  compared 
with  the  result  predicted  by  the  near-wall  diffusion  method. 
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5 . 2 Radiative  Heat  Transfer  Rate 


Radiative  heat  flux  is  also  predicted  by  two  different 
approaches  in  this  study.  From  the  diffusion  approximation, 
the  radiative  flux  is  estimated  by  the  temperature  gradient 
from  the  wall  to  the  adjacent  cell  and  the  radiative 
conductivity  defined  in  Section  2.2.1.  The  expression  is 
written  as 


where  the  temperature  gradient  at  the  wall  is  numerically 
calculated  in  the  same  manner  as  in  the  previous  section. 
For  the  approximation  when  using  the  one-dimensional 
equation  of  radiative  transfer,  the  radiative  flux  into  the 
wall  is  the  net  radiative  heat  flux  resulting  from  the 
exchange  between  the  fluxes  which  come  from  the  inner  gas 
layer  within  5 mean  free  path  (mfp)  and  the  fluxes  emitted 
from  the  wall.  The  net  radiative  flux  at  the  wall  is 
calculated  by 


(5.5) 


2 


(5.6) 


where  the  subscript  w indicates  the  wall. 


CHAPTER  6 

RESULTS  AND  DISCUSSION 

In  this  chapter  the  computational  results  for  the  fluid 
flow  and  heat  transfer  are  presented.  This  chapter  is 
composed  of  two  major  parts,  the  first  of  which  describes 
the  results  for  the  flow  in  both  a straight  tube  and  a tube 
with  a circular  arc  bump.  The  straight  tube  is  utilized  to 
verify  the  fluid  flow  and  heat  transfer  models  used  in  this 
study.  The  tube  with  a circular  arc  bump  is  used  to  test 
the  model's  capability  of  capturing  shock  waves.  The  second 
part  presents  the  numerical  results  for  converging-diverging 
nozzle  flows.  Air  and  UF4  gas  are  the  fluids  considered  in 
this  study.  Air  is  used  to  compare  the  numerical  results 
with  existing  experimental  or  analytical  data,  and  UF4  is 
used  to  simulate  the  flow  of  a UTVR-MHD  system. 

6 • 1 Calculations  in  a Straight  Tube  and  a Tube 
with  a Circular  Arc  Bump 

6.1.1  Verification  of  Modelling 

The  developed  model  is  validated  in  this  section.  Two 
simple  geometrical  configurations,  a straight  tube  and  a 
tube  with  a circular  arc  bump  on  the  inner  wall,  are  used  to 
verify  the  present  model.  The  straight  tube  has  a diameter- 
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to-length  (L/D)  ratio  of  60,  and  a 98  x 60  grid  is  used. 

The  velocity  distributions  within  fully-developed  regions 
are  presented  for  Reynolds  numbers  within  a range  of  1.6  x 
104  < Re  < 3 . 9 x 105.  These  numerical  results  are  compared 
to  the  velocity  profile  depicted  by  the  l/7th  power  law  in 
Figure  6.1.  The  velocity  distribution  of  the  l/7th  power 
law  represents  the  fully-developed  turbulent  flow 
distribution  [24].  As  shown  in  Figure  6.1,  the  results  are 
in  reasonable  agreement.  The  comparison  is,  however,  not 
suited  for  a wide  range  of  Reynolds  numbers. 

Therefore,  the  universal  velocity  distribution,  which 
represents  the  velocity  distribution  for  fully-developed 
turbulent  tube  flows  in  a wide  range  of  Reynolds  numbers,  is 
used  for  the  verification  of  turbulent  flow  solution.  The 
numerical  results  are  compared  to  Laufer's  and  Deissler's 
experimental  data  [24]  for  the  mean  velocity  in  a fully- 
developed  isothermal  tube  flow  as  shown  in  Figure  6.2.  The 
excellent  agreement  between  the  present  model  and  the 
experimental  data  indicates  the  validity  of  the  developed 
model  and  the  related  assumptions  for  prediction  of 
turbulent  flows. 

In  addition,  a close-up  view  of  the  streamwise 
direction  velocity  profiles  at  various  locations  downstream 
are  presented  in  Figure  6.3.  The  temperature  profiles, 
which  are  obtained  at  a constant  wall  temperature  of  1000  K 
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Figure  6.1  Calculated  velocity  distribution  in  a straight  tube  for 
varying  Reynolds  number  and  velocity  distribution 
following  l/7th  law. 
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Figure  6.2  Velocity  distribution  for  turbulent,  isothermal  flow  in 
tubes . 
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Figure  6.3  Close-up  view  of  a typical  near-wall  velocity  profile 
as  calculated  for  an  isothermal  situation. 


62 


and  an  inlet  stagnation  temperature  of  1500  K,  are  presented 
in  Figure  6.4.  For  both  results,  parabolic  profiles  of 
boundary  layers  are  observed.  In  Figure  6.5,  the  pressure 
drop  along  the  axis  is  shown.  The  pressure  drop  is 
presented  as  the  ratio  of  static  to  stagnation  pressure. 

All  of  the  above  results  are  obtained  from  the  solution 
of  the  flows  at  the  steady  state.  In  order  to  verify  the 
steady  state  solution,  the  convergence  history  for  the 
velocity  residual  is  examined,  as  shown  in  Figure  6.6.  The 
steady  state  is  reached  when  the  velocity  residual 
approaches  an  asymptotic  value.  In  this  test  case,  the 
steady  state  is  reached  in  around  80  iterations  with  a 
residual  of  about  4 x 10~3. 

Another  geometry  of  interest  to  this  study  is  a 
converging-diverging  nozzle  in  which  transonic  flow  and 
shock  waves  occur . The  problem  of  flow  in  a channel  with  an 
inward  circular  arc  bump  is  chosen  to  evaluate  the  code  for 
subsonic  and  transonic  steady-state  modeling.  This 
particular  problem  is  well  suited  for  testing  of  solutions 
of  the  subsonic  and  transonic  flow.  The  geometry  and  the 
grids  are  easy  to  generate  accurately,  and  the  symmetry  and 
geometrical  simplicity  aid  in  the  interpretation  of  the 
results . 

For  subsonic  and  transonic  modeling,  the  geometrical 
configuration  is  shown  in  Figure  6.7,  and  the  thickness  of 
the  bump  is  10%  of  the  radius  of  the  channel.  The 
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Figure  6.4  Close-up  view  of  a typical  near-wall  temperature 

profile  as  calculated  for  an  isothermal  situation. 
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Figure  6.5  Calculated  surface  pressure  distribution  along  the 
tube . 
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Figure  6.6  Convergence  rate  for  the  steady  state  solution  in 
straight  tube. 


66 


Center  Line 


Inlet 


Tube  Wall 


Restriction 


1.0 


0.1  | 


3.0 


Outlet 


1.0 


1.0 


Figure  6.7  Geometry  of  a tube  with  a circular  restriction. 


Figure  6.8  A 90x60  grid  system. 
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computational  mesh,  composed  of  90  x 60  grid  lines  as  shown 
in  Figure  6.8,  is  used  in  both  subsonic  and  transonic 
calculations.  With  the  geometrical  configuration  mentioned 
above,  the  boundary  conditions  for  subsonic  flow  are  as 
follow.  Solid-wall  boundary  and  symmetric  boundary 
conditions  are  applied  to  the  solid  walls  and  the  center 
line.  The  inflow  boundary  is  on  the  left  side.  The 
stagnation  pressure  at  the  inlet  is  1 atm,  while  the  back 
pressure  is  kept  at  0.8  atm.  The  inlet  temperature  is 
maintained  at  2500  K,  and  the  wall  temperature  is  kept  at 
1000  K. 

Figure  6.9  displays  the  isomach  lines,  the  normalized 
pressure,  and  the  normalized  temperature  lines  at  steady 
state.  As  depicted  in  the  Figure  6.9',  the  numerical 
solutions  are  quite  symmetric  about  the  midchord,  which  is  a 
good  indication  of  the  accuracy  of  the  solution  for  this 
subsonic  application.  The  surface  pressure  distribution 
along  the  inner  tube  wall  plotted  in  Figure  6.10  presents 
the  symmetric  characteristic  quite  well. 

For  the  case  of  transonic  flow,  the  same  flow 
conditions  used  in  the  subsonic  flow  are  applied  but  the 
back  pressure  is  changed  to  0.25  atm.  From  the 
calculation,  a supersonic  region  appears  in  the  solution 
which  is  depicted  by  a shock  as  shown  in  Figure  6.11. 

The  transient  isomach  contours,  recorded  at  t=0.0003, 
0.0008,  0.0013,  and  0.0022  seconds,  respectively,  are  shown  in 
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(a) 


(b) 


(c) 


Figure  6.9  Computed  contours  at  steady  state  for  subsonic  flow. 

(Tln=2500  K,  Tw=1000  K,  Pm=l  atm,  Pout=0 . 8 atm,  and 
Pr=0.72)  (a)  Mach  number  (b)  normalized  pressure 

(c)  normalized  temperature 
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Figure  6.10  Computed  pressure  distributions  for  a subsonic  flow 
regime;  surface  pressure  (solid  line)  and  centerline 
pressure  (dotted  line)  (Tln=2500  K,  Tw=1000  K,  Pin=l 
atm,  Pout=0.8  atm,  and  Pr=0.72). 
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Figure  6.11.  The  isomach  contour  of  the  steady-state 
solution  shows  a clear  sharp  shock  wave  which  is  located 
around  the  trailing  edge  of  the  bump. 

The  normalized  pressure  and  temperature  contours  are 
also  presented  in  Figures  6.12  and  6.13,  respectively.  The 
same  shapes  of  shock  formations  are  observed.  In  addition, 
the  clear  discontinuity  due  to  the  shock  waves  is  observed 
in  the  surface  pressure  distribution,  as  shown  in  Figure 
6.14. 

6.1.2  Effect  of  Boundary  Cell  Size  in  Convective  Heat 
Transfer  Rate 

As  previously  mentioned,  convective  heat  transfer  at 
the  wall  occurs  by  molecular  diffusion  and  is  determined  by 
Fourier's  law.  When  using  Fourier's  law  to  determine  the 
convective  wall  heat  flux,  the  size  of  the  boundary  grids  is 
of  crucial  importance  because  variation  in  the  size  of  the 
boundary  grids  significantly  influences  the  calculation  of 
the  temperature  gradient  which  is  directly  proportional  to 
the  wall  heat  flux. 

A mesh  refinement  study  is  conducted  to  determine  the 
largest  grid  size  below  which  the  computed  wall  heat 
transfer  rate  does  not  introduce  a significant  error. 

Figure  6.15  presents  the  convective  heat  transfer  rates 
estimated  as  a function  of  the  size  of  the  boundary  grid, 
dy,  in  a 60  x 60  grid  system.  The  dy  is  the  distance  from 
the  wall  to  the  point  where  the  radial  temperature  gradient 
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Figure 


(b) 


(d) 


.11  Calculated  isomach  number  lines  at  different  time  for 
transonic  flow.  (Tln=2500  K,  Tw=1000  K,  Pin=l  atm, 
Pout=0.25  atm,  and  Pr=0.72)  (a)  t=0.0003  sec.  (b) 
t=0 .0008  sec.  (c)  t=0.0013  sec.  (d)  t=0.0022  sec. 
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Figure  6.12  Calculated  normalized  pressure  contours 
time  for  transonic  flow.  (Tln=2500  K,  Tw= 
atm,  Pout=0.25  atm,  and  Pr=0.72)  (a)  t=0. 

t=0.0008  sec.  (c)  t=0.0013  sec.  (d)  t=0 


at  different 
=1000  K,  Pin— 1 
0003  sec.  (b) 
.0022  sec. 
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(d) 


Figure  6.13  Calculated  normalized  temperature  contours  at 
different  time  for  transonic  flow.  (Tln=2500  K,  T, 
K,  Pm=l  atm,  Pout=0.25  atm,  and  Pr=0.72)  (a)  t=( 

sec.  (b)  t=0.0008  sec.  (c)  t=0.0013  sec.  (d)  t=( 
sec . 


=1000 

.0003 

.0022 
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Figure  6.14  Calculated  pressure  distributions  for  a transonic  flow 
regime;  surface  pressure  (solid  line)  and  centerline 
pressure  (dotted  line)  (Tln=2500  K,  Tw=1000  K,  Pm=l 
atm,  Pout=0.25  atm,  and  Pr=0.72). 
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Figure  6.15  Heat  transfer  rates  obtained  by  Navier-Stokes  solver 
for  various  boundary  cell  sizes.  A 60x60  grid  is  used. 
(Tln=4000  K,  Tw=1800  K,  Pm=l  atm,  and  Pout=0 . 5 atm) 
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at  the  wall  is  calculated.  More  detail  of  the  grid  sizes 
near  the  wall  is  given  in  Table  6.1.  In  this  calculation,  a 
straight  tube,  5 cm  in  diameter  and  1 m in  length,  is  used. 
The  initial  inlet  temperature  is  4000  K,  and  the  wall 
temperature  is  kept  at  1800  K.  The  stagnation  pressure  at 
the  inlet  is  1 atm,  while  the  back  pressure  is  0.5  atm. 

For  a grid  size,  dy,  greater  than  2.91  x 10“5  m,  the 
computed  heat  flux  is  underestimated,  as  is  shown  in  Figure 
6.15.  For  the  grid  sizes  of  dy  = 2.67  x 10-6  and  dy  = 1.03 
x 10"5,  the  calculated  heat  transfer  rates  are  equal.  In 
terms  of  the  dimensionless  distance  from  the  wall,  y+  = 
yv*/v,  where  v*  is  the  frictional  velocity  and  v is  the 
kinematic  velocity,  the  maximum  allowable  boundary  grid  size 
is  approximately  5.  This  indicates  that  accurate  prediction 
of  the  temperature  gradient  near  the  wall  requires  grid 
sizes  less  than  the  thickness  of  the  laminar  sublayer.  Note 
that  the  dimensionless  thickness,  y+,  of  the  laminar 
sublayer  in  turbulent  flow  is  5 [40]  . In  this  case,  as 
shown  in  Table  6.1,  seven  node  points  are  located  within  the 
laminar  sublayer,  y+  = 5. 

A grid  sensitivity  analysis  is  performed  to  confirm  the 
relation  between  the  thickness  of  the  laminar  sublayer  and 
the  maximum  allowable  grid  size  near  the  wall.  The 
thickness  of  the  laminar  sublayer  is  varied  with  the  use  of 
different  flow  conditions,  such  as  system  pressures  or  inlet 
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Table  6.1  Details  of  grid  sizes  near  wall  for  a tube  in  5 cm 

diameter  and  1 m length.  60  x 60  grid  system.  Tin=4000 
K,  Twaii=1800  K,  Pin=l  atm,  and  Pout=  0.5  atm. 


Grid  # 
from  Wall 

Actual  Distance  from  Wall 
in  m (x  10+5) 

y 

<^> 

V 

1 

0.266 

0.83 

2 

0.693 

1.32 

3 

1.035 

1.88 

4 

1.421 

2.52 

5 

1.858 

3.24 

6 

2.353 

4.06 

7 

2.914 

4.98 

8 

3.548 

6.02 

9 

4.266 

7.20 

10 

5.079 

8.54 

11 

5.998 

10.05 

12 

7.040 

11.77 
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stagnation  temperatures.  For  instance,  higher  inlet 
stagnation  pressure  and  back  pressure  are  chosen  to 
calculate  the  convective  wall  heat  flux  while  other 
conditions  are  kept  the  same. 

Figure  6.16  shows  the  convective  wall  heat  transfer 
rates  which  are  estimated  for  different  choices  of  grid 
size.  The  inlet  and  back  pressures  in  this  case  are  10  and 
9.5  atm,  respectively,  and  the  grid  size  is  60  x 80. 

However,  the  inlet  and  exit  temperatures  are  kept  the  same 
as  those  in  the  60  x 60  grid  case. 

In  this  case  the  estimated  heat  fluxes  are  more 
sensitive  to  the  grid  size,  as  shown  in  Figure  6.16. 
Therefore,  it  is  evident  from  Figure  6.16  that  if  the  first 
grid  from  the  wall  is  located  in  the  laminar  sublayer,  the 
estimated  convective  heat  transfer  rate  is  accurate  and 
reasonable.  Table  6.2  lists  the  calculated  grid  sizes  of  a 
60  x 80  grid  system.  According  to  the  data  given  in  Table 
6.2,  two  grid  points  are  located  within  the  laminar 
sublayer,  even  though  more  grids  in  the  radial  direction  are 
allocated.  This  indicates  that  the  thickness  of  the  viscous 
sublayer  is  thinner  for  the  high  pressure  case  where  the 
higher  density  increases  the  value  of  Reynolds  numbers. 

From  this  analysis,  it  is  evident  that  accurate 
convective  heat  transfer  rates  can  be  estimated  if  the 
boundary  cell  size  in  the  radial  direction  is  smaller  than 
the  laminar  sublayer,  y*<5,  thickness.  The  use  of  a boundary 
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Table  6.2  Details  of  grid  sizes  near  wall  for  a tube  in  5 cm 

diameter  and  1 m length.  60  x 80  grid  system.  Tln=4000 
K,  Twall=1800  K,  Pln=10  atm,  and  Pout=9.5. 


Grid  # 
from  Wall 

Actual  Distance  from  Wall 
in  m (x  10+5) 

y+ 

(#) 

V 

1 

0.125 

1.65 

2 

0.386 

3.46 

3 

0.672 

5.43 

4 

0.983 

7.58 

5 

1.323 

9.92 

6 

1.694 

12.48 

7 

2.098 

15.28 

8 

2.540 

18.32 

9 

3.022 

21.65 

10 

3.548 

25.28 

11 

4.122 

29.25 

12 

4.749 

33.57 
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Figure  6.16  Heat  transfer  rates  obtained  by  Navier-Stokes  solver 
for  various  boundary  cell  sizes.  A 60x80  grid  is  used. 
(Tln=4000  K,  Tw=1800  K,  Pm=10  atm,  and  POUL=9.5  atm) 
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cell  size  larger  than  the  laminar  sublayer  thickness  results 
in  an  underestimate  of  the  calculated  wall  heat  flux. 

Using  the  present  heat  transfer  model,  the  calculated 
convective  wall  heat  transfer  rate  is  compared  with  existing 
experimental  data  for  a developing,  isothermal  pipe  flow. 

The  result  is  presented  in  Figure  6.17  and  compared 
with  Boelter' s experimental  data  [12]  for  the  air  flow  at  a 
Reynolds  number  of  53000  at  the  exit.  The  agreement  between 
the  calculated  result  and  the  experimental  data  is 
reasonably  good.  Moreover,  these  results  agree  better 
beyond  the  diameter-to-length  ratio  (L/D)  of  about  10. 

Figure  6.18  presents  the  results  of  the  comparison 
between  the  directly  calculated  Nusselt  number  from  the 
transverse  temperature  gradient  of  the  thermal  field 
solution  and  the  Nusselt  number  evaluated  using  the  Dittus- 
Boelter  correlation  [2] . For  this  test  calculation,  the 
stagnation  inlet  and  the  back  pressure  are  10  atm  and  9.5 
atm,  respectively,  and  the  inlet  temperature  is  4000  K and 
the  wall  temperature  is  maintained  at  1800  K.  Results  are 
in  reasonably  good  agreement  as  the  flow  becomes  fully 
developed.  From  this  analysis,  it  is  concluded  that  the 
convective  wall  heat  transfer  rate  is  accurately  predicted 
by  the  developed  fluid  flow  and  heat  transfer  model. 


82 


CM  CM 


O 

o 

CM 


o 

id 


o 

o 


o 

id 


o 

o 


;ui(nN)/nN 


Q 

\ 

INI 


Figure  6.17  Nusselt  number  versus  axial  position  for  a developing 
isothermal  pipe  flow  at  a Reynolds  number  of  53000. 
(Nu)  lnf  is  the  Nusselt  number  evaluated  at  Z/D=20. 
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Figure  6.18  Comparison  of  Nusselt  number  as  obtained  from  the 

computed  heat  flux  with  the  result  obtained  by  using 
the  Dittus-Boelter  correlation.  A 60x60  grid  is  used. 
(Tln=4 000  K,  Tw=l 800  K,  Pm=10  atm,  and  Pout=9 . 5 atm) 
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6.1.3  Evaluation  of  Diffusion  Approximation 

An  assessment  is  made  of  the  diffusion  approximation  as 
a matter  of  interest.  The  local  radiative  heat  flux  at 
various  tube  wall  locations  is  computed  by  two  different 
approximations,  the  diffusion  and  the  simplified  integral 
approximation.  A comparison  is  made  using  the  same  tube 
dimension  as  that  used  in  the  previous  section  with  a 60  x 
60  grid.  The  inlet  temperature  is  taken  to  be  4000  K while 
the  tube  wall  temperature  is  assumed  to  be  1800  K. 

As  previously  mentioned,  the  Rosseland's  diffusion 
approximation,  which  can  be  applicable  for  a certain  range 
of  absorptivities,  is  used  to  account  for  the  radiative  heat 
transfer.  The  limiting  range  of  the  absorptivity  for  the 
diffusion  model  is  examined  by  comparing  the  radiative  heat 
flux  estimated  by  the  two  different  approximations  mentioned 
above . 

The  comparative  results  are  plotted  in  Figure  6.19 
which  presents  the  ratio  of  the  heat  flux  of  the  diffusion 
model  to  that  of  the  integral  model  at  various  axial 
positions.  This  comparison  is  performed  by  varying  the 
absorptivity  with  a change  of  flow  field  pressure.  The  mean 
absorptivities,  based  on  a tube  wall  temperature  of  1800  K 
and  inlet  pressures  ranging  from  10  atm  to  100  atm,  vary 
from  a value  of  111.1  cm"1  to  a value  of  1111.1  cm"1.  As 
seen  in  Figure  6.19,  the  diffusion  approximation 
consistently  underestimates  the  heat  flux  by  about  30% 
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Figure  6.19  Comparative  result  for  radiative  heat  flux  between 

diffusion  approximation  and  1-D  integral  approximation 
for  varying  gas  opacity  due  to  different  flow 
conditions . 
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compared  to  the  results  of  the  integral  approximation  for 
mean  absorptivities  greater  than  555.5  cm'1.  Moreover,  it 
is  found  that  the  discrepancy  between  the  diffusion 
approximation  and  the  simplified  integral  approximation 
becomes  more  significant  when  applied  to  an  optically  thin 
gas,  such  as  a gas  with  an  absorptivity  of  less  than  about 
100  cm'1.  Therefore,  it  is  suggested  that  the  diffusion 
approximation  should  not  be  used  for  an  optically  thin  gas. 

6.1.4  Convective  and  Radiative  Heat  Transfer 

The  heat  fluxes  obtained  with  the  developed  models  are 
presented  for  both  convective  and  radiative  heat  transfer. 
The  convective  and  radiative  heat  transfer  rates  presented 
in  this  section  are  calculated  by  using  the  near-wall 
diffusion  method  and  the  1-D  integral  approximation, 
respectively.  The  prediction  of  the  radiative  heat  transfer 
including  the  effect  of  the  tube  wall  temperature  and  the 
gas  inlet  temperature  are  presented  here.  In  order  to 
examine  the  wall  temperature  effect,  the  tube  wall 
temperatures  of  choice  are  1800,  2100,  and  2400  K,  while  the 
inlet  temperature  is  kept  at  4000  K. 

The  heat  fluxes  obtained  with  the  different  wall 
temperatures  are  plotted  in  Figure  6.20.  It  is  seen  that 
the  radiative  fluxes  increase  as  higher  wall  temperatures 
are  selected.  This  behavior  is  due  simply  to  the  fact  that 
the  radiative  wall  heat  flux  is  dependent  on  the  local 
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Figure  6.20  Radiative  heat  flux  calculated  for  varying  wall 

temperatures.  (Tln=4000  K,  Pin=10  atm,  and  PouC— 9.5  atm) 
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conditions  at  the  wall  as  well  as  the  gas  inlet 
temperatures.  This  result  is  displayed  in  Figure  6.21.  If 
the  gas  inlet  temperatures  are  high,  the  high  radiative 
fluxes  occur  at  the  wall,  as  expected,  but  the  effect  of 
changing  the  inlet  gas  temperature  is  less  than  that  of 
changing  the  wall  temperature.  This  indicates  that  the 
radiative  heat  transfer  at  the  wall  is  more  likely  to  be 
affected  by  the  local  conditions  at  the  wall. 

Using  the  developed  model,  the  heat  fluxes  due  to  the 
combined  radiative  and  convective  heat  transfer  for  the  five 
cases  summarized  in  Table  6.3  are  presented.  Figures  6.22 


Table  6.3  Flow  conditions  for  calculation  and  comparison  of 
radiative  and  convective  heat  transfer  rates. 


Pin  (atm) 

Pout  (atm) 

Tin  ( K ) 

T«  ( K ) 

CASE  1 

10 

9.5 

4000 

1800 

CASE  2 

10 

9.5 

3000 

1800 

CASE  3 

10 

9.5 

4000 

2100 

CASE  4 

30 

29.5 

4000 

1800 

CASE  5 

50 

49.5 

4000 

1800 

and  6.23  display  the  results  of  the  calculations  for  cases  1 
and  2 that  differ  only  in  the  gas  inlet  temperature  which 
also  changes  the  value  of  the  mean  opacity  of  the  gas.  The 
radiative  flux  in  Figure  6.22,  for  an  inlet  temperature  of 
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Figure  6.22  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (Tln=4000  K,  Tw=1800  K,  Pin=10  atm,  and 
Pout=9.5  atm) 
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Figure  6.23  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (T,n=3000  K,  Tw=1800  K,  Pm=10  atm,  and 
Pout  = 9.5  atm) 
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4000  K,  is  almost  two  times  larger  than  the  convective  flux, 
while  the  radiative  flux  for  the  inlet  temperature  of  3000  K 
is  slightly  less  than  two  times  in  the  magnitude  of  the 
convective  flux.  Both  radiative  and  convective  fluxes  for 
case  1 are  almost  twice  as  large  as  those  of  case  2 where 
the  inlet  gas  temperature  is  reduced  by  1000  K. 

Results  of  radiative  and  convective  heat  transfer 
calculations  for  case  3 are  presented  in  Figure  6.24.  Flow 
conditions  in  case  3 are  similar  to  those  of  case  1,  with 
the  exception  that  the  wall  temperature  is  raised  by  300  K. 
It  is  seen  that  the  radiative  flux  is  about  60%  greater  than 
that  of  case  1,  while  the  convective  flux  remains  almost  the 
same . 

Finally,  with  the  change  of  the  flow  field  pressure, 
the  radiative  and  convective  heat  fluxes  are  estimated  and 
plotted  in  Figures  6.25  and  6.26.  The  flow  conditions  for 
cases  4 and  5 are  similar  to  those  of  case  1,  but  the  flow 
field  pressure  is  increased  by  factors  of  3 and  5.  As  can 
be  seen  in  the  figures  mentioned  above,  the  convective  heat 
flux  increases,  but  the  radiative  heat  flux  decreases  as  the 
flow  field  pressure  becomes  higher  and  higher.  The  increase 
of  the  convective  heat  flux  seems  to  follow  the  change  of 
the  density  which  is  due  to  the  pressure  increase.  The 
decrease  in  the  magnitude  of  the  radiative  flux  may  be 
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Figure  6.24  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (Tln=4000  K,  Tw=2100  K,  Pm=10  atm,  and 
Pout=9.5  atm) 
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Figure  6.25  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (Tln=4000  K,  Tu=1800  K,  Pm=30  atm,  and 
Pou,=29.5  atm) 
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Figure  6.26  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (Tln=4000  K,  Tw=1800  K,  Pin=50  atm,  and 
Pout=4  9.5  atm) 
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simply  caused  by  the  shortening  of  the  radiation  mean  free 
path  by  the  higher  molecular  density.  It  is  interesting  to 
note  that  the  magnitudes  of  the  convective  heat  flux  for 
cases  4 and  5 are  greater  than  those  of  the  radiative  flux. 

6.1.5  Effect  of  Internal  Heat  Generation  in  Heat  Transfer 
Rate 

UTVR-MHD  systems  are  characterized  by  high  internal 
heat  generation  rates  in  the  gas.  The  developed  model  is 
used  to  predict  the  heat  transfer  rates  for  the  flow  of 
gases  with  internal  heat  generation.  Uniform  power 
densities  ranging  from  0 to  1000  w/cc  are  assumed  and 
included  in  the  energy  equation  as  source  terms. 

The  inlet  gas  temperature  for  the  case  under 
consideration  is  2000  K which  is  a typical  temperature  at 
the  core  inlet  of  a UTVR-MHD  system.  The  wall  temperature 
is  kept  at  1800  K,  and  the  inlet  stagnation  pressure  is  10 
atm.  The  back  pressure  is  maintained  at  9.5  atm. 

First,  the  evolution  of  bulk  temperature  as  a function 
of  the  tube  length  for  different  values  of  power  density  is 
shown  in  Figure  6.27.  The  increase  of  the  bulk  temperature 
at  the  power  density  of  1000  w/cc  is  about  2000  K.  This 
indicates  that  a power  density  higher  than  1000  w/cc  is 
required  to  maintain  the  exit  gas  temperature  of  the  tube  at 
about  4000  K.  The  temperature  of  4000  K is  a typical 
temperature  at  the  core  exit  required  to  achieve  high 
efficiency  in  the  MHD  energy  conversion  system. 
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Figure  6.27  Bulk  temperature  evolution  for  power  densities  of  0, 
50,  100,  500  w/cc,  and  1000  w/cc.  (Tln=2000  K,  Tw=1800 
K,  P in—  1 0 atm,  and  PouC=9.5  atm) 
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Figure  6.28  presents  the  Nusselt  numbers  for  the 
convective  heat  transfer  at  different  power  densities.  The 
Nusselt  numbers  for  power  densities  of  50  w/cc  and  100  w/cc 
progressively  increase  in  regions  beyond  the  tube  entrance. 
However,  for  high  power  densities  such  as  500  w/cc  and  1000 
w/cc,  a drastic  increase  of  Nusselt  numbers  occurs  near  the 
tube  entrance. 

These  results  indicate  that  for  moderate  power 
densities,  50  w/cc  and  100  w/cc,  the  heat  generated  due  to 
the  internal  heat  generation  and  the  heat  removed  by 
convection  heat  transfer  in  the  axial  direction  are  almost 
comparable  for  the  positions,  L/D  <5.  As  the  flow  proceeds 
downstream,  the  generated  heat  steadily  exceeds  the  heat 
transferred  through  the  wall.  For  very  high  power  densities 
such  as  500  w/cc  and  1000  w/cc,  the  heat  from  the  internal 
heat  source  considerably  exceeds  the  transferred  heat 
through  the  wall  as  the  flow  proceeds  only  a short  distance 
into  the  tube. 

Radiative  heat  transfer  rates  are  estimated  for  power 
densities  of  10  w/cc,  100  w/cc,  and  500  w/cc,  and  the 
results  are  plotted  in  Figure  6.29.  Since  the  radiative 
heat  transfer  is  more  strongly  dependent  on  the  gas 
temperature,  the  radiative  heat  transfer  rate  increases 
rapidly  along  the  tube  at  a power  density  of  500  w/cc. 
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6.1.6  Test  of  Convergence 

Since  an  iterative  procedure  is  used  to  obtain  a 
steady-state  solution  as  a time  asymptotic  solution  of  the 
unsteady  equations,  a test  of  convergence  is  necessary. 

Thus,  the  convergence  of  the  process  is  examined  for 
different  flow  conditions.  The  velocity  residuals  as  a 
function  of  the  iteration  numbers  are  plotted  in  Figure 
6.30,  and  a steady  state  is  reached  when  the  velocity 
residual  approaches  an  asymptotic  value.  As  seen  in  these 
figures,  the  number  of  iterations  for  the  steady  state 
solutions  is  increased  as  pressure  and  temperature  increase. 

6 . 2 Calculations  in  a Converging-Diverging  Nozzle 
6.2.1  Solution  of  Flow  and  Thermal  Field  in  a Nozzle 

The  axisymmetric  Navier-Stokes  code  with  the  Baldwin- 
Lomax  turbulent  model  is  used  to  carry  out  a simulation  of 
the  converging-diverging  nozzle  in  which  the  flow  field  is 
under  the  influence  of  large  pressure  gradients.  With  the 
thermal  field  solution  obtained,  convective  and  radiative 
heat  transfer  rates  along  the  wall  of  the  nozzle  are 
predicted.  Since  the  simplified  integral  approach,  which  is 
applicable  only  to  a simple  tube  geometry,  is  inaccurate  for 
the  prediction  of  the  radiative  heat  transfer  for  this 
nozzle  geometry,  the  diffusion  approximation  is  adopted. 


Velocity  Residual  Velocity  residual 
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(a) 


(b) 


Figure  6.30  Convergence  rates  for  different  flow  conditions. 

(a)  Tln=4000  K,  Tw=1800  K,  Pm=10  atm,  and  Pouc=9.5  atm, 
Q=0 

(b)  Tln=5000  K,  Tw=1800  K,  Pm=10  atm,  and  Pout=9.5  atm, 
Q=0 


Velocity  Residual  Velocity  Residual 
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(d) 


Figure  6.30 — continued. 

(c)  Tln=4000  K,  Tw=1800  K,  Pm=30  atm,  and  Pout=29.5  atm, 
Q=0 

Tln=4000  K,  Tw=1800  K,  Pin=10  atm,  and  Pout=9.5  atm, 
Q=50  MW 
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The  nozzle  configuration  used  is  shown  in  a schematic 
diagram  in  Figure  6.31,  and  the  flow  conditions  utilized  are 
as  follows:  The  inlet  stagnation  pressure  is  20  atm;  the 
inlet  gas  temperature  is  4000  K;  the  constant  wall 
temperature  is  kept  at  1800  K;  and  UF„  gas  properties 
evaluated  at  2000  K are  used.  A mesh  is  generated  that  has 
60  points  in  the  streamwise  direction  and  60  points  in  the 
radial  direction  for  the  lower  half  of  the  nozzle  due  to  the 
geometric  symmetry  about  the  center  line.  In  the  radial 
direction  the  mesh  is  compressed  near  the  solid  surface  of 
the  lower  boundary  to  resolve  the  viscous  layer.  The  upper 
boundary  of  the  simulation  is  the  nozzle  centerline.  Figure 
6.32  shows  the  grid  system  mentioned  above.  Results  of  the 
simulation  are  presented  in  terms  of  contours,  velocity 
vector  plots,  pressure  distribution  and  heat  transfer  rates 
along  the  solid  wall. 

Figure  6.33  depicts  the  flow  pattern  which  indicates 
that  the  flow  is  accelerated  through  the  throat  of  the 
nozzle.  A close  view  of  the  velocity  and  the  temperature  at 
various  locations  along  the  solid  wall  are  presented  in 
Figures  6.34  and  6.35.  The  hydrodynamic  and  thermal 
boundary  layers  are  observed  in  both  figures.  Figures  6.36 
and  6.37  show  the  Mach  number  contours  and  the  pressure 
contours  for  the  nozzle.  The  subsonic  inlet  flow 
accelerates  from  left  to  right  and  reaches  sonic  speeds  near 


1 


0.36  m 


Figure  6.31  A converging-diverging  nozzle 
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Figure  6.32  A 60x60  grid  in  a converging-diverging  nozzle. 
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Figure  6.33  Flow  pattern  in  a converging-diverging  nozzle 
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Figure  6.34 


Close-up  view  of  velocity  distribution. 
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1 2 345  673 


Figure  6.35  Close-up  view  of  temperature  distribution. 
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(a) 


(b) 


Figure  6.36  contours. 

(a)  isomach  number  contour. 

(b)  normalized  pressure  contour. 
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the  throat  and  continues  to  accelerates  as  the  flow  proceeds 
downstream.  As  can  be  seen  from  both  the  pressure  and  the 
Mach  number  contour,  compression  waves  are  initiated  along 
the  wall  surface  just  downstream  from  the  throat  toward  the 
axis  of  symmetry.  These  waves  coalesce  to  form  a weak 
oblique  shock  which  is  further  propagated  as  it  approaches 
the  axis.  Moreover,  the  presence  of  an  oblique  shock  at  an 
axial  location  just  downstream  from  the  throat  is  clearly 
shown  in  the  normalized  pressure  distribution  in  which  the 
static-to-stagnation  pressure  ratio  is  plotted  as  a function 
of  axial  distance,  L/Dthrc,  along  the  nozzle  in  Figure  6.37. 

6.2.2  Convective  and  Radiative  Heat  Transfer 

The  variation  of  the  heat  fluxes  along  the  nozzle  are 
shown  in  Figure  6.38  for  the  same  flow  conditions  given  in 
the  previous  section  and  the  nozzle  shown  in  Figure  6.31. 
Figure  6.38  reveals  that  the  maximum  value  of  the  convective 
heat  flux  occurs  just  upstream  of  the  throat  in  the  vicinity 
of  the  mass  flux  maximum.  The  radiative  heat  flux  exhibits 
a trend  similar  to  that  of  the  convective  heat  flux  because 
both  fluxes  are  proportional  to  the  same  temperature 
gradient  at  the  wall.  With  the  selection  of  higher 
stagnation  pressures,  50  atm  and  100  atm,  for  the  inlet 
boundary  condition,  it  is  evident  from  Figure  6.39  that  the 
heat  fluxes  increase,  as  expected,  with  increasing 
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Figure  6.37  Normalized  pressure  distribution  along  the  nozzle 
surface . 
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Figure  6.38  Convective  and  radiative  heat  fluxes  for  flow 

conditions.  (Tln=4000  K,  Tw=1800  K,  Pm=20  atm,  and 
Pout=l  atm) 
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Figure  6.39  Convective  heat  flux  for  different  stagnation 
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stagnation  pressures  which  correspond  to  larger  mass  fluxes. 
The  variation  of  heat  fluxes  along  the  nozzle  as  a function 
of  changes  in  geometrical  parameters  of  the  nozzle,  such  as 
a throat  curvature  and  an  uphill  angle  of  the  nozzle,  is 
also  examined.  The  obtained  heat  fluxes  decrease  with  an 
increase  of  the  uphill  angle  toward  the  entrance  of  the 
nozzle,  as  shown  in  Figure  6.40,  but  there  is  no  significant 
change  of  the  maximum  heat  flux  at  the  throat.  Figure  6.41 
presents  results  for  changes  in  the  throat  curvature.  As 
shown  in  Figure  6.41,  the  heat  fluxes  are  only  broadened  by 
increasing  the  throat  curvature,  while  the  peak  value  of  the 
heat  flux  at  the  throat  remains  unchanged.  Thus,  it  can  be 
stated  that  changes  in  the  geometrical  parameters  do  not 
affect  the  maximum  heat  flux  at  the  throat.  In  all  the 
figures  of  the  heat  fluxes,  small  peaks  in  the  middle  of  the 
upstream  region  are  observed.  These  peaks  are  a result  of 
the  existence  of  weak  secondary  flows. 

6.2.3  Effect  of  Grid  System 

As  previously  mentioned,  a disadvantage  of  algebraic 
methods  for  grid  generation  is  that  they  tend  to  yield 
skewed  grid  systems,  which  may  affect  the  flow  and  thermal 
field  solution  and,  furthermore,  may  create  a more  severe 
problem,  such  as  a convergence  problem  in  a time-marching 
scheme.  Thus,  the  effect  of  the  grid  system  for  the 
solution  at  the  steady  state  is  investigated  by  comparing 
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Figure  6.41  Convective  heat  flux  for  different  radii  of  curvature 
at  throat . 
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results  obtained  from  both  the  skewed  and  orthogonal  grid 
systems.  The  skewed  and  orthogonal  grid  system  for  a bell- 
type  nozzle  are  displayed  in  Figures  6.42  and  6.43.  Flow 
conditions  are  the  same  as  those  of  the  previous  section, 
but  a different  choice  of  grid  size,  40x60,  is  used. 

Figures  6.44  (a)  and  (b)  represent  the  flow  patterns  for  the 
skewed  and  the  orthogonal  grid  system,  respectively.  An 
additional  comparison  is  performed  with  the  contours  of 
pressure,  density,  velocity,  and  Mach  number.  The  contours 
for  both  grid  systems  are  presented  in  Figures  6.45  through 
6.48.  From  these  figures,  it  is  observed  that  the  results 
obtained  from  the  different  grid  systems  are  almost  the 
same.  Further,  the  static-to-stagnation  pressure  ratios 
along  the  nozzle  in  Figures  6.49  (a)  and  (b)  are  equal.  In 
conclusion,  there  is  no  advantage  in  using  an  orthogonal 
grid  system  over  using  a skewed  grid  system  for  the  simple 
nozzle  configuration  examined  in  this  work. 

6.2.4  Assessment  for  Using  Frequency  Averaged  Rosseland 
Absorption  Coefficient  in  Radiative  Heat  Transfer 

In  this  study,  frequency  averaged  absorption 
coefficients  are  used  to  predict  radiative  transfer.  Using 
the  mean  absorption  coefficient  eliminates  carrying  out  a 
spectral  analysis  and  integrating  over  all  wavelengths  to 
obtain  total  energy  quantities.  The  question  is  whether  the 
analysis  will  yield  an  accurate  solution  for  a particular 
problem.  Thus,  an  analysis  is  performed  to  assess  the 
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Figure  6.42  A skewed  grid  in  a bell-type  nozzle. 
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Figure  6.43 


An  orthogonal  grid  in  a bell-type  nozzle. 
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(b) 


Figure  6.44  Flow  patterns. 

(a)  skewed  grid 

(b)  orthogonal  grid 
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(a) 


Figure  6.45  Normalized  pressure  contours. 

(a)  skewed  grid 

(b)  orthogonal  grid 
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(a) 


(b) 


Figure  6.46  Normalized  density  contours. 

(a)  skewed  grid 

(b)  orthogonal  grid 


(b) 


Figure  6.47  Normalized  velocity  contours. 

(a)  skewed  grid 

(b)  orthogonal  grid 
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(b) 


Figure  6.48  Mach  number  contours. 

(a)  skewed  grid 

(b)  orthogonal  grid 
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Figure  6.49  Surface  normalized  pressure  distribution  along  the 
nozzle . 

(a)  skewed  grid 

(b)  orthogonal  grid 
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sensitivity  of  radiative  transfer  to  the  mean  absorption 
coefficient.  The  absorption  spectrum  of  UF6  from  2000  to 
4200  A [31]  is  used.  The  analysis  is  performed  with  a mean 
absorption  coefficient  over  all  the  spectral  range  and  mean 
absorption  coefficients  for  two  wavelength  bands,  2000  A- 
3000  A and  3000  A-4200  A.  The  Rosseland  mean  opacity  over 
all  wavelengths,  2000  A-4200  A,  is  5.0  x 103  (1/m)  for 
temperature,  4000  K,  and  pressure,  10  atm.  The  Rosseland 
mean  opacities  for  the  wavelength  bands  of  2000  A-3000  A 
and  3000  A-4200  A are  1.06  x 104  (1/m)  and  3.8  x 10  (1/m), 
respectively . 

Using  the  above  opacities,  radiative  heat  fluxes  are 
estimated  with  the  diffusion  approximation  and  1-D  integral 
approximation.  As  expected,  the  diffusion  approximation  is 
very  sensitive  to  the  magnitude  of  the  opacity.  For  the 
wavelength  range  of  2000  A-3000  A,  the  radiative  heat  flux 
is  underpredicted  as  plotted  in  Figure  6.50.  For  the  range 
of  3000  A-4200  A,  the  diffusion  approximation  is  not 
applicable  because  the  criterion  for  the  diffusion 
approximation,  described  in  subsection  2.2.1,  is  not 
satisfied  by  the  long  mean  free  path  of  photon.  As  seen  in 
this  analysis,  the  Rosseland  mean  absorption  coefficient 
over  all  wavelengths  may  have  a large  value,  but  the 
spectral  absorption  coefficient  may  be  very  small  in  certain 
spectral  regions.  Therefore,  it  is  recommended  to  use 
wavelength  bands  in  which  the  spectral  absorption 
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Figure  6.50  Radiative  heat  transfer  rates  calculated  by  the 

diffusion  approximation  for  a mean  opacity  over  all 
wavelengths  and  for  a mean  opacity  over  the  range  o 
2000  A-3000  A. 
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coefficient  is  everywhere  large  and  to  evaluate  a Rosseland 
mean  for  each  of  these  regions. 

The  1-D  integral  approximation  is  less  sensitive  to  the 
magnitude  of  the  Rosseland  mean  opacity.  Radiative  heat 
transfer  rates  are  calculated  with  the  mean  opacity  over  all 
wavelengths  and  the  mean  opacities  for  the  ranges  of  2000  A- 
3000  A and  3000  A-4200  A.  The  result  for  the  mean  opacity 
is  compared  with  the  averaged  heat  flux  for  the  opacities  of 
two  wavelength  bands.  Figure  6.51  shows  the  results  and 
reveals  that  use  of  the  mean  opacity  for  1-D  integral 
approximation  can  be  acceptable  for  the  prediction  of 
radiative  heat  transfer. 
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CHAPTER  7 

SUMMARY  AND  CONCLUSIONS 


An  investigation  into  the  use  of  computational  fluid 
dynamics  (CFD)  has  been  carried  out  to  predict  the  flow 
field  solutions  and  expected  heat  transfer  rates  of  interest 
to  Ultrahigh  Temperature  Vapor  Core  Reactor  systems  using  a 
straight  tube  geometry  and  the  converging-diverging  geometry 
of  the  reactor  exit  nozzle.  An  axisymmetric  two-dimensional 
(2-D)  Navier-Stokes  based  model,  which  is  composed  of  an 
algebraic  two-layer  eddy  viscosity  turbulence  model,  a 
convective  heat  transfer  model,  and  radiative  heat  transfer 
models,  is  developed.  Using  the  developed  models,  the 
following  studies  have  been  accomplished. 

1.  The  flow  and  thermal  fields  at  a high  Reynolds 
number  and  at  a high  temperature  are  expressed  in  a thin- 
layer  Navier-Stokes  equations.  The  governing  equations  are 
solved  by  a hybrid  implicit-explicit  method  based  on  the 
finite  volume  approach.  The  numerical  algorithm  used  in 
this  study  has  a rapid  convergence  compared  to  a simple 
explicit  method.  Even  though  very  fine  grids  are  used  for 
the  prediction  of  convective  heat  flux  into  the  wall,  less 
than  300  iterations  are  required  to  reach  the  steady-state 
for  any  flow  conditions  taken.  The  computed  velocity 
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distribution  in  a straight  tube  for  air  is  compared  to 
existing  experimental  data.  Good  agreement  between  the 
computed  results  and  the  measured  data  indicates  the 
validity  of  the  developed  model  and  the  related  assumptions 
for  prediction  of  flows. 

2.  Applicability  for  transonic  flow  modeling  is 
examined  by  checking  the  capability  of  capturing  shock  waves 
in  a channel  with  an  inward  circular  arc  bump.  The  computed 
data  of  the  Mach  number  and  pressure  are  presented  in  Mach 
number  contours  and  normalized  pressure  contours,  which  show 
clear,  sharp  shock  waves  located  around  the  trailing  edge  of 
the  bump.  Applicability  for  transonic  flow  modeling  is  also 
assured  by  observing  the  discontinuity  due  to  a shock  wave 
in  the  surface  pressure  distribution. 

3.  Convective  heat  flux  into  the  wall  is  directly 
determined  by  the  known  temperature  gradients  at  the  wall  if 
extremely  fine  grids  contain  a grid  point  within  the  laminar 
sublayer,  y+  - 5.  The  computed  heat  transfer  rates  in  the 
Nusselt  number  are  compared  to  the  measured  data  for  air  and 
the  result  obtained  with  an  empirical  correlation,  Dittus- 
Boelter  correlation.  Good  agreements  between  the  computed 
result  and  the  measured  data  and  the  result  of  the  empirical 
correlation  indicate  that  the  developed  flow  and  heat 
transfer  model  is  accurate. 

4.  Radiative  heat  flux  is  predicted  by  both  the 
diffusion  approximation  and  the  one-dimensional  equation  of 
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radiative  transfer.  From  the  comparison  of  the  two  results 
obtained,  it  is  found  that  the  radiative  heat  flux  obtained 
by  the  diffusion  approximation  is  underestimated  by  about 
30%  compared  to  the  other  results  for  mean  absorptivities 
greater  than  555.5  cm-1,  evaluated  at  the  stagnation 
pressure,  50  atm,  and  the  wall  temperature,  1800  K. 

5.  The  effect  of  internal  heat  generation  on  heat 
transfer  in  straight  tubes  is  examined  for  a variety  of 
power  densities,  0 to  1000  w/cc.  For  very  high  power 
densities  such  as  500  w/cc  and  1000  w/cc,  the  heat  supplied 
by  the  internal  heat  sources  considerably  exceeds  the  heat 
transferred  through  the  wall. 

6.  The  convective  and  radiative  heat  fluxes  are 
predicted  for  a converging-diverging  nozzle.  The  maximum 
value  of  convective  and  radiative  heat  fluxes,  as  expected, 
occurs  just  upstream  of  the  throat  in  the  vicinity  where  the 
mass  flux  is  a maximum.  The  heat  fluxes,  in  general, 
increase  with  increasing  stagnation  pressures  which 
correspond  to  larger  mass  fluxes. 

The  above  accomplishments  demonstrate  that  the 
computational  fluid  dynamics  and  heat  transfer  model  based 
on  the  axisymmetric,  thin-layer,  Navier-Stokes  equations 
represents  a viable  computational  tool  for  predicting  flow 
and  thermal  fields  of  interest  pertaining  to  UTVR  systems. 
Moreover,  the  model  developed  can  be  applied  to  many 
engineering  applications  related  to  energy  system,  such  as 
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combustion  chambers,  Gas  Core  Reactors  (GCR) , and  nuclear 
propulsion  systems. 

Convective  heat  flux  is  accurately  predicted  by  the 
convective  heat  transfer  model,  but  an  assessment  of 
radiative  heat  transfer  models  in  comparison  with 
experimental  data  is  not  made  because  experimental  data  of 
radiative  heat  transfer  rates  for  an  optically  thick  gas  are 
not  available.  Therefore,  it  is  desirable  that  hardware  be 
built  to  test  the  radiative  heat  transfer  rates  obtained  by 
the  radiative  heat  transfer  models. 


APPENDIX 

SOURCE  LISTING  OF  CODE 


This  is  the  source  program  of  a Navier-Stokes  based  model  for 
the  analysis  of  the  combined  radiative  and  convective  heat 
transfer. 


c 

c 


program  main 


common/ a/rho (100,100,2), rhou (100,100,2) , rhov (100, 100,2) , 
c e (100, 100,2) , fsl (2) , fs2 (2) , 

c u(100,  100) , v ( 100 , 100)  ,ei (100,  100) , p ( 1 00 , 100) , 

c fs3 (2) ,fs4 (2) , z (101, 101)  , r (101, 101) , vol (101, 101) , 

c gamma, gml , ggml, cv, pr, vou,  ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm,n,m, kl, k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

common/b/iturb, prt, ep (100, 100) 
common/c/du (100,100,4), dus (100,100,4) 
common/e/gll (101),gl2(101),gl3(101),gl4(101), 

1 g21 (101) ,g22 (101) ,g23 (101) ,g24 (101) , 

2 g31 (101) ,g32 (101) ,g33 (101) ,g34 (101) , 

3 g41  (101)  ,g4 2(101)  ,g4 3(101)  ,g44  (101)  , 

4 ff (101) , implt 
common /spit/ if lxs, sz, sr 
common/q/qp 

common /gm/thrt 

dimension  ynu (100) , znum (100) , rr (100) 
common/rad/nradq, rabsco 

dimension  pu  (100, 100), pv (100, 100), newtw (100) 
dimension  rnr (100), ul 7(100), ul 10(100) 
common/tem/t (101,101), ttwall 
common/speed/rmach (100, 100) 




c gamma  - ratio  of  specific  heats,  cp/cv 
c cv  - specific  heat  at  constant  volume 
c pr  - prandtl  number 
c prt  - turbulent  prandtl  number 
c iflxs  = 0 no  flux  splitting 
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c =1  first  order  flux  splitting  included 

c =2  second  order  flux  splitting  included 

c il  - range  of  mesh  index  i 

c jl  - range  of  mesh  index  j 

c nend  - range  of  time  step  index  nadv 

c nprnt  - print  parameter  (complete  flow  field  printed  at 
c nadv=nend, surface  properties  printed  at  nadv=k*nprnt , 

c for  k=l, 2, 3 . . . , nprnt  must  be  a divisor  of  nend) 

c cfl  - courant-f riedrich-lewy  number  (must  be  less  than  one) 

c iread  = 0 - new  start 

c = 1 - start  and  continue  from  stored  file 

c iwrite  = 0 - no  data  file  to  be  written  for  restarting 
c = 1 - data  file  to  be  written  for  restarting 

c implt  = 0 - explicit  maccormack  method 
c > 0 - implicit  gauss  seidel  line  relaxation 

c iturb  = 0 - laminar  flow,  = 1 - turbulent  flow 
c nvisc  - value  of  nadv  at  which  viscous  flow  computation  begins 
c jlfm  - estimate  of  range  of  j index  containing  boundary  layer 

c iadbwl  = 1 - adiabatic  wall,  = 0 - isothermal  wall 

c twall  - temperature  of  wall 

c ptt  - total  pressure  specified  at  entrance 

c ttt  - total  temperature  specified  at  entrance 

c vou  - tangent  of  flow  angle  at  entrance 

c pe  pressure  specified  at  exit 

c qp  - internal  heat  source (w/m**3) 
c nradq  = 0 - no  radiation  heat  transfer 
c = 1 - radiation  heat  transfer 

c rabsco  - rosseland  absorption  coefficient 
c igrid  = 0 - skew  grid 
c = 1 - orthogonal  grid 

c 

c 

open  (unit=15, file='  tube . inp' , status=' old'  ) 
open (unit=16, f ile='  tube . out' , status=' old' ) 
open (unit=22, f ile=' univ . dat ' , status=' old' ) 
open (unit =50, file='  graph.dat' , status=' old' ) 
open (unit=999, f ile=' subl . dat' , status=' old' ) 
open (unit =7 7, f ile='  heat 2 . dat' , status=' old' ) 
open  (unit =60,  f ile='  resdu . dat'  , status='  old'  ) 
c 

c gamma  is  1.4  for  air  and  gamma  is  1.08  for  UF4. 
c gamma=1.4 

gamma=l . 08 
gml=gamma-l . 0 
ggml=gamma*gml 

c cv  is  319.31  for  UF 4 property  estimated  at  2000  k in  j/kg  k. 
c cv  is  717.5  for  air  in  j/kg  K. 

cv=31 9 . 31 
c cv=717 . 5 

c Prandtl  number  of  uf4  is  0.915. 

c Prandtl  number  of  air  is  0.72 
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c 


c 


c 

c 

c 


c 

c 

c 


pr=0. 915 
pr=0 . 72 
prt=0 . 90 

read  (15,*)  il, jl, nend, nprnt, iread,  iwrite, iflxs, igrid 
read  (15,*)  implt , iturb, nvisc, jlfm, nradq, rabsco, ttwall 
read  (15,*)  iadbwl, twall, cf 1, vou, ptt, ttt , pe, qp 

il=il+mod (il, 2) 
ilml=il-l 
jlml= jl-1 
eiwall=cv* twall 
do  1 i=l, il 
do  1 j=l , jl 
ep (i, j) =0 . 0 
1 continue 

Initial  imformations  are  stored  for  graphics. 

write (50,*)  nend 
write (50,*)  nprnt 
write (50, *) il 
write (50, *) jl 
write (50, *) m 
write (50, *) ptt 

if  (iread. eq.0)  call  gmtry 

call  prntxy 

time=0 . 0 
nadv=0 

if  (iread. eq.0)  call  initl 

if  (iread. eq.l)  call  rdwrt (1, time) 

nstart=nadv+l 

lcnt=l 

write  (16,5) 

5 format (lhl) 
write  (16,10) 

10  format (///, 56x, 18hi  n p u t data,//) 
write  (16,20)  il, iturb, cfl 

20  format (lh  ,30x,8hil  =, 2x, i8 , 7x, 8hiturb  =,2x,i8,7x, 

c 8hcf 1 =, 2x, f 8 . 3) 

write  (16,30)  jl, iadbwl, if lxs 
30  format (lh  ,30x,8hjl  =, 2x, i8, 7x, 8hiadbwl  =,2x,i8,7x, 

c 8hiflxs  =,2x,i8) 
write  (16,40)  nend, twall, nstart 
40  format (lh  , 30x, 8hnend  =, 2x, i8, 7x, 8htwall  =, 2x, f 8 . 3, 7x, 
c 8hnstart  =,2x,i8) 
write  (16,50)  nprnt , nvisc, implt 
50  format (lh  ,30x,8hnprnt  =, 2x, i8 , 7x, 8hnvisc  =,2x,i8,7x, 
c 8h  implt  =, 2x, i8) 
write  (16,80) 
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80  format (////, 51x, 28hupstream  boundary  conditions) 
write  (16, 90) 

90  format  (/,  32x,  19x,  3hptt,  19x,  3httt,  19x,  3hvou) 
write  (16,100)  ptt,ttt,vou 
100  format (32x, 3 (lOx, fl2. 4) ) 
m=l 
n=2 
c 

call  be 

if  (nend. eq. 0 . and. iwrite . eq. 1)  call  rdwrt (2, time) 
if  (nend.eq.0)  stop 
kl=l 
k2=2 
c 

do  111  i = 1, il 
do  111  j = 1, jl 
pu (i, j) =0 . 0 
pv(i,  j)=0.0 
111  continue 
c 

do  200  nadv=nstart , nend 
do  140  i=l , il 
do  140  j=l, jl 

si=sqrt ( (z (i, j+1 ) -z (i, j ) ) **2  + (r (i, j+1) -r (i, j) ) **2) 
s zp=+ (r (i, j + 1) -r (i, j) ) /si 
srp=- ( z (i, j+1) -z (i, j) ) /si 
uu=sxp*u  (i, j) +syp*v (i, j) 

sj=sqrt ( (z (i+1, j)-z (i, j) ) **2+(r (i+1, j)-r (i, j) ) **2) 

sxp=-(r(i  + l,  j)-r(i,  j))/sj 

syp=+ (z (i  + 1, j) -z  (i, j) ) /s j 

vv=szp*u (i, j) +srp*v (i, j) 

cc=sqrt (ggml *abs (ei (i,  j)  ) ) 

dti=0 . 5*vol (i, j) /si/ (abs (uu) +cc) 

dt j=0 . 5*vol (i, j) /s j/ (abs (vv) +cc) 

if  (i . eq. 1 . and. j . eq. 1 ) then 

dte=dti 

dt=cf l*dti 

endif 

if  (dti.lt.dte)  dte=dti 
if  (dtj.lt. dte)  dte=dt j 
140  continue 

if  (implt.eq.0)  dt=cfl*dte 

if (mod (nadv, nprnt ) .eq.0)  write  (16, 150)  nadv,dte,dt 
150  format (lOx, 5hnadv=, i7, 5x, 4hdte=, el2 . 4, 5x,  3hdt=,  el2 .4) 
if  (nadv. ge.nvisc. and. iturb. eq. 1)  call  turbfl 
time=time+dt 
call  1 (dt ) 
sum=0 . 0 

do  300  i =1 , il 
do  300  j =1 , j 1 

sum=sum+sqrt ( (u(i, j)-pU(i,  j) ) **2  + (v (i, j) -pv (i, j) ) **2) 

continue 


300 
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sum=sum/ (il* jl) 
do  400  i =1, il 
do  400  j =1 , jl 
pu  (i,  j)  =u  (i,  j) 
pv(i,  j)  =v  (i,  j) 

400  continue 

c if (mod (nadv, nprnt ) . eq. 0 . or . nadv . eq. nend) then 

do  155  i=l,il 
do  155  j=l , jl 
t (i,  j)  =ei (i,  j) / cv 
155  continue 

c 
c 

c - calculate  mach  number 
do  157  i=l , il 
do  157  j=l , jl 
c=sqrt (ggml*abs (ei (i,  j)  ) ) 

rmach  (i,  j)  =sqrt  (u  (i,  j)  *u  (i,  j)  +v  (i,  j)  *v  (i,  j)  ) /c 
157  continue 

if  (mod (nadv, nprnt) .eq. 0)  write (16, 160)  nadv, dt, time 
write ( 6, * ) ' iter  #=' , nadv, ' , residual=' , sum 

write (60,*)  nadv, sum 

160  format (lOx,  5hnadv=, i7, 5x, 3hdt=,  el2 . 4, 5x,  5htime=, el2 . 4) 

if  (mod(nadv, nprnt) .eq.0. or. nadv. eq. nend)  call  prntff(time) 
200  continue 


c calculation  of  heat  flux 


first  order  error  heat  flux  calculation 


dd=thrt *2 . 
ilm3=il-3 
do  153  i=2,il-l 
dr=0 . 5* (r (i,  3)-r (i,  2)  ) 
dt=t (i, 2) -twall 
dtt=t (i, 2) -twall 
dtdr=dt /dr 
dttdr=dtt /dr 

drl=(r(i,  4)  — r (i,  2)  ) +0 . 5*  (r  (i,  5) -r  (i,  4)  ) 
dr2= (r (i,  8) -r (i, 2) ) +0.5* (r (i, 9)-r (i, 8) ) 
dr3= (r (i, 12) -r (i, 2) ) +0 . 5* (r (i,  13) -r (i, 12) ) 
dl=dr 

d2=  (r(i,3)-r(i,2)  ) + 0.5*  (r(i,4)-r(i,3)  ) 
d3=drl 

d4=  (r(i,5)-r(i,2)  ) +0.5*  (r(i,  6)  -r  (i,  5)  ) 

d5=(r  (i,  6)-r  (i,2)  ) +0.5*  (r  (i,  7)-r  (i,  6)  ) 

d6=(r(i,7)-r(i,2)  ) +0.5*  (r(i,8)-r(i,7)  ) 
d7=dr2 

d8=(r  (i,  9)-r  (i,2)  ) +0.5*  (r  (i,  10) -r  (i,  9)  ) 
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d9= (r  (i, 10) -r (i, 2) ) +0.5* (r(i, 11) -r  (i, 10) ) 
dl0= (r (i, 11) -r  (i, 2) ) +0.5* (r(i, 12) -r (i, 11) ) 
dll=dr3 

dl2=  (r  (i,  13)  -r  (i,  2)  ) +0 . 5*  (r  (i,  14) -r (i,  13)  ) 
d20=(r  (i,20)-r  (i,2)  )+0.5*(r(i,21)-r(i,20)) 
dtl=t (i, 4 ) -twall 
dt2=t (i, 8) -twall 
dt3=t (i, 12) -twall 
dt20=t (i, 20) -twall 
dtdrl=dtl/drl 
dtdr2=dt2/dr2 
dtdr3=dt3/dr3 
dtdr20=dt20/d20 
do  1110  j=l , jl-1 
if ( j .eq. 1)  then 
rr ( j) =0 .5* (r (i, 3) -r (i, 2) ) 
else 

rr(j)  = (r  (i,  j+l)-r(i,2)  )+0.5*<r  (i,  j+2)-r(i,  j + 1)  ) 
endif 

1110  continue 

if (i .eq. ilm3) then 

write (6, *) ' 1, dr' , i, dr, drl, dr2, dr 3 
write (6, *) ' dl  to  d6' , dl, d2, d3, d4, d5, d6 
write (6, *) ' d7  to  dl2' , d7,  d8,  d9,  dlO,  dll,  dl2 

endif 


c - air  viscosity  (Sutherland's  formula) 

c rmu=l . 458e-06  * sqrt  (t  (i,  2)  **3)  / (t  (i,  2) +110 . 4) 

c - uf4  viscosity  at  2000  K 

c rmu=3 . 357e-6*sqrt (t (i,  2) ) / (aa+bb*t (i, 2) ) 
rmu=8 . 67e-5 

c - thermal  conductivity  rk=cv*gamma*rmu/pr 

rk=cv*gamma*rmu/pr 
zod=0 . 5* (z (i,2) +z (i+1,2)  ) /dd 


c second  order  error  heat  flux  calculation 


c dd=2 . *thrt 

c do  153  i=2, il-2 

c drmins=0.5* (y (i, 3) -y (i, 2) ) 

c dr=0.5*  (r  (i,4)-r(i,3)  ) 

c drplus=dr+drmins 

c alpa=drplus/drmins 

c - calculate  radiative  heat  flux 
if (nradq.eq. 1) then 
sigma=5 . 67e-8 

rabsco-rosseland  absorption  coefficient- 
tt=0 . 5* (t (i, 2) +t  (i, 3) ) 

radk= (8 . *sigma* (ttwall**3) ) / (3 . *rabsco*2 . ) 
endif 


c 

c 


141 


c dtdy=  (-alpa*  (alpa+2. ) *twall+  ( (alpa+1  .)**2)*t(i,2)-t(i,3))/ 
c 1 (alpa* (alpa+1 . ) *dymins ) 

c 

c - qd  is  the  variable  of  convective  heat  transfer  rate  at  the 
c wall. 

c - qrad  is  the  variable  of  radiative  heat  transfer  rates  at  the 
c wall 

c 

qd=l . e-6*rk*dtdy 
if (nradq.eq. 1) then 

qrad=l . e-6*radk*dtdy 
qtot=qd+qrad 

zod=0 .5*  (x(i,2)+x(i+l,2)  ) /dd 
write (77,*)  zod, qd, qrad, qtot 
endif 

153  continue 

c - if  1-D  integaral  approxiamation  is  used  to  calculate  the 
c radiative  heat  transfer  rate,  this  parts  must  be  included, 

c aaw=rabsco* (ttt/ttwall) 

c do  2222  i =2,  il-1 

c plusin=0.0 

c do  2223  j =2,  jl-2 

c rw=r (i, j) -r (i, 2) 

c if (rw. le . (5 . /aaw) ) then 

c rr ( j) =0 . 5* (r (i, j+2) -r  (i, j) 

c ri=( (t (i, j) +t (i, j+1) ) *0.5) **4 

c plusw (i, j) =rabsco*sigma*  (1 . /pi) *ri*exp (-rabsco*rr ( j) ) *rr ( j) 

c plusin=plusin+plusw (i, j) 

c endif 

c2223  continue 

c zod=0.5* (z (i,2)+z (i+1,2) ) /dd 

c plusi (i) =plusin 

c radw (i) =1 .e-6* (plus (i) - (1 . /pi) * sigma* (twall) **4) ) 

c when  1-D  integral  approximation  is  used,  make  an 

c open  statement 

c write (1000,  *)  zod,  radw(i) 

c2222  continue 
c 

if  (iwrite.eq. 1)  call  rdwrt (2, time) 

stop 

end 




subroutine  gmtry 



c 

c reads  computational  mesh  coordinates 
c 

common/ a/ rho (100,100,2) , rhou (100, 100, 2) , rhov (100, 100, 2) , 
c e(100, 100,2)  ,fsl  (2)  ,fs2(2)  , 

c u (100, 100) , v (100, 100) ,ei (100, 100) , p (100, 100) , 
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c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml , ggml, cv,  pr, vou, ptt , ttt , pe, 

c il, jl, ilml, jlml , jlfm, n,  m,  kl,  k2, 

c iadbwl, eiwall,  lent , cfl, 

c nadv, nend, nvisc 

common /gm/thrt 
common /grid/igr id 
dimension  rr(100),  o(100) 
c 

pi=3. 14159 

c - mks  unit  is  used, 
c 

rl=l . 0 
thrt=0 . 025 
c 

c - For  a converging-diverging  nozzle,  rad  is  a radius  of  c 
c curvature  at  throat, 

c 

c rad=0.02 

c - theta  is  uphill  angle  and  epsln  is  downhill  angle, 
c 

c theta=20 . 0*pi/180 . 0 

c 

c - For  a straight  tube,  both  angles  are  zero. 
theta=0 . 0 

c - epsln=20 . 0*pi/180 
epsln=  0.0 
c 

ithrt=il/2+l 
sint=sin (theta) 
cost=cos (theta) 
sine=sin  (epsln) 
cose=cos (epsln) 
dr=thrt /l . Oe+4 
c 

deta=l . 0/  ( jl-2) 

call  stretch (thrt, dy, deta, ck, ekml) 
c 

c cd=ck*deta 

do  1 j=2, jl 
c 

r (ithrt, j) =-thrt* (1.0- (exp ( ( j-2) *ck*deta) -1.0) /ekml) 
c 

1 continue 
c 
c 

dz=l . 0*rl/  ( i 1— 2 ) 
detap=l .0/  ( ithrt-2 ) 

call  stretch (0 . 5*rl, dz,  detap, ckp, ekmlp) 
do  2 i=2 , il 
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if  (i.le.ithrt)  then 

z (i, 2) =0 . 5*rl* (1.0- (exp ( (ithrt-i) *ckp*detap) -1.0) /ekmlp) 
else 

z (i, 2) =0 . 5*rl* (1.0+ (exp ( (i-ithrt) *ckp* detap) -1.0) /ekmlp) 
endif 

if  (z  (i, 2)  . le . 0 . 5*rl-rad*sint ) 
c yri,  2)  =-thrt-rad+rad*cost 

c - (0 . 5*rl-rad*sint-z (i,  2) ) *sint/cost 

if  (z  (i,  2)  . gt . 0 . 5*rl-rad*sint . and. z (i,  2)  . It . 0 . 5*rl+rad*sine) 
c r (i, 2)  = (- (thrt+rad) +sqrt (rad* rad- (z (i,2)-0.5*rl) **2) ) 
if  (z  (i, 2)  . ge . 0 . 5*rl+rad*sine) 
c r (i,  2)  =-thrt-rad+rad*cose 

c - (z (i, 2) -0 . 5*rl-rad*sine) * sine /cose 

do  2 j=3, jl 
z (i,  j)  =z  (i,  2) 

r (i,  j) =-r (i, 2) *r (ithrt, j) /thrt 
2 continue 
c 

c set  boundary  cells  by  linear  extrapolation 
c 

do  10  i=2,il 

z (i,  l)=2.0*z  (i,  2)  -z  (i,3) 
r(i,l)=2.0*r(i,2)-r(i,3) 
z (i, jl  + l)=2.0*z  (i, jl) -z (i, j 1 — 1 ) 
r(i, jl+1 ) =2 . 0*r (i, jl) -r  (i, j 1-1 ) 

10  continue 
c 
c 

jlpl=jl+l 
do  20  j = l, jlpl 

z(l,  j)=2.0*z(2,  j)-z(3,  j) 
r (1,  j)  =2 . 0*r  (2,  j)-r  (3,  j) 
z (il+1, j)=2.0*z  (il, j ) — z (il-1, j) 
r (il+1,  j)=2.0*r(il,  j)  -r  (il-1,  j) 

20  continue 
c 

c calculate  orthogonal  grid  system  by  an  algebraic  technique, 
if  (igrid. eq. 0)  then 
go  to  41 
else 

do  300  i =3,  ithrt-1 

teta=atan ( (r (i  + 1, 2) -r  (i, 2) ) / (z (i  + 1, 2) -z (i, 2) ) ) 
rr  (i)  = (r  (i,  jl)  -r  (i,  2)  ) / sin  (teta) 
o(i)=z  (i,  jl)  + (r(i,  jl)  -r  (i,  2)  ) / tan  (teta) 
do  300  j =3, jl 

phi=atan( (r (i+1, j) -r (i, j) ) / (z (i+1, j) -z (i, j) ) ) 
z (i,  j) =o (i) -rr (i) *cos  (phi) 
r(i,j)=-(o(i)-z(i,j))  *tan  (phi) 

300  continue 

do  301  i=ithrt+l,  il-1 

teta=atan( (r (i  + 1, 2) -r (i, 2) ) / (z (i  + 1, 2) -z  (i, 2) ) ) 
rr(i)  = (r(i,  jl)  — r (i,  2)  ) /sin  (teta) 
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o(i)=z  (i,  jl)  + (r  (i,  jl)  -r  (i,  2)  ) /tan (t eta) 
do  301  j =3, jl 

phi=atan ( (r  (i+1,  j)  -r  (i,  j)  ) / (z  (i+1,  j)  -z  (i,  j)  ) ) 
z (i,  j)=o(i)-rr(i)  *cos  (phi) 
r(i,  j)  =-  (o(i)-z  (i,  j)  ) *tan  (phi) 

301  continue 

endif 

41  continue 

c 

c calculate  all  volumes 

c 

do  30  j=l , jl 
do  30  i=l , il 

vol (i, j)  = (abs ( (z (i,  j)-z (i  + 1, j) ) *r (i+1, j + 1) 

1 + (z  (i  + 1,  j) -z  (i  + 1,  j + 1)  ) *r  (i,  j) 

2 + (z (i+1,  j + 1) -z (i, j) ) *r  (i  + 1, j) ) 

3 +abs ( (z (i,  j ) - z (i  + 1, j + 1) ) *r(i, j + 1) 

4 + (z  (i+1,  j + 1) -z  (i,  j + 1)  ) *r  (i,  j) 

5 +(z (i, j+l)-z (i, j) ) *r (i+1, j+1) ) ) /2.0 
30  continue 

c 

c reset  the  volumes  of  the  boundary  cells  equal  to  that  of 
c their  flowfield  counterparts, 
c 
c 

do  40  i=l,il 

vol (i, 1) =vol  (i, 2) 

vol (i, jl) =vol (i, jlml) 

40  continue 
return 
end 
c 



subroutine  stretch (rr,  dd,  deta, ck, ekml) 




ck=l . 0 

do  10  k=l, 10 
c do  10  k=l, 18 

ek=exp (ck) 
ekdeta=exp (ck*deta) 
fk=dd-rr* (ekdeta-1 . 0) / (ek-1 .0) 

fkp-~rr* (deta*ekdeta- (ekdeta-1 .0) *ek/ (ek-1 . 0) ) / (ek-1 . 0) 
ck=ck-fk/fkp 
10  continue 

ekml=exp (ck) -1 . 0 

return 

end 
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subroutine  initl 



c 

c form  initial  flow  field 

c 

common/a/rho (100, 100,2), rhou (100, 100,2), rhov (100, 100,2), 
c e (100,  100, 2) , fsl  (2) , fs2  (2) , 

c u (100,  100) , v (100, 100) ,ei(100,100) , p (100,  100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101, 101) , 

c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il,  jl,  ilml, jlml, jlfm, n, m, kl, k2 , 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

common/ e/gll (101) ,gl2 (101),gl3(101),gl4(101), 
c g21 (101) ,g22 (101) ,g23 (101) ,g24 (101)  , 

c g31 (101) ,g32 (101) ,g33(101)  ,g34 (101)  , 

c g41 (101) ,g42 (101) ,g43 (101)  ,g44 (101)  , 

c ff (101) , implt 

c 

eit=cv*ttt 
rtt=ptt/ (gml*eit) 
do  10  j=l , jl 
do  10  i=l,il 
rho (i, j, 1) =rtt 
rhou (i, j, 1) =0 . 0 
rhov (i, j, 1) =0 . 0 
e (i, j, 1) =rtt*eit 
if  (i. eq.il)  then 
rho (i, j,  1) =rtt*pe/ptt 
e (i, j, 1) =pe/gml 
endif 

u (i,  j)  =rhou  (i,  j,  1)  /rho  (i,  j,  1) 
v (i,  j)_=rhov  (i,  j,  1)  /rho  (i,  j,  1) 

ei  (i,  j)  =e  (i,  j,  1)  / rho  (i,  j,  1)  -0 . 5*  (u  (i,  j)  **2+v  (i,  j)  **2) 

p (i, j) =gml*rho (i, j, 1) *ei  (i, j) 

rho  (i,  j,  2)  =rho  (i,  j,  1) 

rhou ( i , j , 2 ) =rhou ( i , j , 1 ) 

rhov ( i , j , 2 ) =rhov ( i , j , 1 ) 

e (i,  j,  2)  =e  (i,  j,  1) 

10  continue 
c 

return 

end 


subroutine  lbc(dt) 


common/ a/ rho (100,100,2), rhou (100,100,2),  rhov (100,100,2), 
c e (100,  100,  2) , fsl  (2) , fs2  (2) , 

c u(100, 100) , v ( 1 00 , 100) ,ei (100, 100) , p ( 1 00 , 100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol  (101, 101) , 
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c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il,jl,ilml,jlml,jlfm,n,m,kl,k2, 

c iadbwl, eiwall, lent , cfl, 

c nadv, nend, nvisc 

common/c/du (100,100,4)  ,dus (100,100,4) 
common/e/gll  (101) ,gl2 (101),gl3(101),gl4(101), 
c g21  (101) ,g22 (101) ,g23 (101) ,g24 (101) , 

c g31  (101) ,g32 (101)  ,g33 (101) , g34 (101)  , 

c g41 (101) ,g42 (101)  ,g43 (101) ,g44 (101)  , 

c ff (101) , implt 

gpl=gamma+l . 0 
c 

c inlet  be  -constant  total  pressure  and  temperature 
c 

do  10  j=2, jlml 
rtt=ptt/ (gml*cv*ttt) 
astr2= (2 . 0/gpl) *ggml*cv*ttt 
rr=rho (1, j, n) 
cc=sqrt (ggml*ei (1, j) ) 
uu=u (1, j) 

dd4=dt* (uu-cc) / (0.5* (z (3, j)— z (1, j) ) ) 

r4=-  (dd4/  (1 . 0-dd4)  ) * (p(2,  j)  -p  (1,  j) -rr*cc*  (u(2,  j)  -u  (1,  j)  ) ) 
hh=ptt* (2 . 0*gamma/gpl) * (uu/astr2) * 
c (1.0- (gml/gpl) *uu*uu/astr2 ) ** (1 . 0/gml) 

um=uu-r4/ (hh+rr*cc) 

rm=rtt* (1.0- (gml/gpl) *um*um/astr2 ) ** (1 . 0/gml) 
eim=cv*ttt*  (1.0- (gml/gpl) *um*um/astr2) 
du  (1,  j,  1) =rm-rho (1, j, n) 
du  (1,  j,  2) =rm*um-rhou (1, j, n) 
du  (1,  j,  3)  =0 . 0 

du ( 1 , j , 4 ) =rm* (eim+0 . 5*um*um) -e ( 1 , j , n) 

10  continue 
c 

c exit  be  subsonic-pe  specified  or  supersonic 
c 

do  20  j=2, jlml 

si=sqrt ( (z  (il, j + l)-z (il, j) ) **2+ (r (il, j + l)-r (il, j) ) **2) 

szp=  (r (il, j+1) -r (il, j) ) /si 

srp=- ( z (il, j+1) -z (il, j) ) /si 

rr=rho (ilml, j, n) 

up=  szp*u (ilml, j) +srp*v(ilml,  j) 

vp=-srp*u (ilml, j) +szp*v (ilml, j) 

cc=sqrt (ggml*ei (ilml, j) ) 

rccsq=l . 0/ (cc*cc) 

dd=si*dt/vol (ilml, j) 

ddl=up*dd 

dd2= (up+cc) *dd 

dd4= (up-cc) *dd 

rl=~ (ddl/ (1 . 0+ddl) ) * (rho (il, j, n) -rho (ilml, j, n) 
c - (p (il, j) -p (ilml, j) ) *rccsq) 

r2=- (dd2/ (1 . 0+dd2) ) * (p (il,  j) -p (ilml, j) 
c +rr*cc* (sxp*u (il, j) +syp*v (il, j) -up) ) 
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r3=- (ddl/ (1 . O+ddl) ) * (-syp*u (il,j)+sxp*v(il,j) -vp) 
r4=- (dd4/ (1 . 0+dd4) ) * (p (il, j ) -p (ilml, j) 
c -rr*cc*(sxp*u(il,j) +syp*v (il, j) -up) ) 

dp=0 .5* (r2+r4 ) 
if  (up.lt.cc)  dp=0.0 
dup= (r2-dp) / (rr*cc) 
dvp=r3 

rm=rho (il, j, n) +rl+dp*rccsq 
um=u (il, j) +szp*dup-srp*dvp 
vm=v (il, j ) +srp*dup+szp*dvp 
du (il,  j,  1) =rm-rho (il, j, n) 
du (il,  j,  2) =rm*um-rhou (il, j, n) 
du (il, j, 3) =rm*vm-rhov (il, j, n) 

du (il, j, 4)  = (p (il, j) +dp) /gml  + O . 5*rm* (um*um+vm*vm) -e  (il, j, n) 
20  continue 
return 
end 
c 



subroutine  be 



c 

common/a/rho (100,  100,2), rhou (100,  100,2),  rhov(100,  100,2), 
c e (100,  100,2) , fsl (2) , fs2 (2) , 

c u(100,  100) , v (100,  100) ,ei (100,  100)  ,p (100,  100) , 

c fs3 (2) , fs4 (2) , z (101, 101)  , r (101,  101)  , vol (101, 101) , 

c gamma, gml, ggml, cv, pr, vou, ptt, ttt, pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl,  k2, 

c iadbwl, eiwall, lent,  cfl, 

c nadv, nend, nvisc 

common/ e/gll (101) ,gl2 (101) ,gl3(101)  ,gl4 (101) , 
c g21 (101) ,g22 (101) ,g23 (101) ,g24 (101) , 

c g31 (101) ,g32 (101) ,g33 (101) ,g34 (101) , 

c g41 (101) , g42 (101), g4 3(101),  g4  4 (101), 

c ff (101) , implt 


gpl=gamma+l . 0 

c **  be  for  inlet  ** 

do  10  j=2, jlml 
rtt=ptt / (cv*gml*ttt) 
astr2= (2 . 0/gpl) *ggml*cv*ttt 
um=u ( 1 , j ) 
v (1, j) =0 . 0 

rho (1 , j , m) =rtt * (1.0- (gml/gpl ) *um*um/astr2 ) ** ( 1 . 0/gml ) 
ei  (1^  j) =cv*ttt* (1.0- (gml/gpl) *um*um/astr2) 

P (1, j) =gml*rho (1, j,m) *ei (1, j) 
rhou ( 1 , j , m) =rho ( 1 , j , m) *um 
rhov (1, j , m) =0 . 0 

e (1, j,m) =rho (1,  j,m) * (ei (1, j ) +0 . 5*um*um) 

10  continue 
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c **  be  for  exit  ** 

pe=p (il, jl) 
if lag=0 

do  20  jz=2, jlml 
j= jlml+2- jz 

si=sqrt ( (z (il, j + 1) -z (il,  j)  ) **2  + (r (il,  j + 1) -r (il,  j)  ) **2) 

szp=  (r (il, j+1) -r (il, j) ) /si 

srp=- (z (il, j+1) -z (il, j) ) /si 

up=szp*u (ilml, j) +srp*v (ilml, j) 

upil=szp*u (il, j)+srp*v(il,  j) 

cc=sqrt (ggml*ei (ilml,  j)  ) 

ccil=sqrt (ggml*ei (il,  j)  ) 

if  (up.gt .cc. and. upil.gt. ceil)  then 

if lag=l 

go  to  20 

endif 

if  (iflag.eq.l)  go  to  15 
dr= (pe-p  (il, j) ) / (cc*cc) 

P (il, j) =pe 

rho (il, j,m) =rho (il, j,m) +dr 
rhou (il, j,m) =rho (il, j,m) *u (il,  j) 
rhov (il, j,m) =rho (il, j,m) *v (il, j) 
ei(il,j)=p(il,j)  / (gml*rho  (il,  j,m)  ) 

e (il,  j,m)  =rho  (il,  j,m)*(ei(il,j)+0.5*(u(il,j)  **2+v  (il,  j) **2) ) 
go  to  20 
15  continue 

rho  (il,  j,m)  =rho  (ilml,  j,m) 
rhou (il, j, m) =rhou (ilml, j,  m) 
rhov ( il , j , m) =rhov ( ilml , j , m) 
e (il,  j,m)  =e  (ilml,  j,m) 
u (il,  j) =u (ilml, j) 
v (il,  j)  =v (ilml, j) 
p (il,  j) =p (ilml, j) 
ei (il,  j)  =ei (ilml, j) 

20  continue 
c 

c lower  wall  and  centerline 
c 

do  30  i=l,il 
u (i,  1) =-u (i, 2) 
v (i,  1) =-v  (i, 2) 
ei (i,  1) =ei (i, 2) 

if (iadbwl . eq. 0)  ei (i,  1) =2 . 0*eiwall-ei  (i, 2) 
p(i,l)=p(i,2) 

rho (i,  1, m) =rho (i, 2,m) *ei (i, 2) /ei  (i, 1) 

u (i, jl) =u (i, jlml) 

v (i, jl) =-v (i, jlml) 

ei (i, jl) =ei (i, jlml) 

p (i,  jl) =p (i, jlml) 

rho (i, jl, m) =rho (i, jlml,m) 

30  continue 
c 


return 

end 
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subroutine  fsi(i,j,ii) 

c 

c 

common/a /rho (100,100,2), rhou (100,100,2),  rhov (100,  100,2) , 
c e (100,  100, 2) , fsl  (2) , fs2  (2) , 

c u (100,  100) , v (100,  100) ,ei (100,  100) ,p (100,  100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol  (101, 101) , 

c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl,  k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

common /b/ it urb, prt, ep (100, 100) 
common/e/gll (101) ,gl2 (101),gl3(101),gl4(101), 
c g21 (101) , g22 (101), g23 (101) , g2 4(101), 

c g31 (101) , g32 (101) , g33 (101) , g34 (101) , 

c g41 (101) , g42 (101) ,g43(101) ,g44 (101) , 

c ff (101) , implt 

common/ spit /if lxs,  sz , sr 
common/rad/nradq,  rabsco 
c 

aa=0 . 80 
bb=-7 . 10e-5 

sz=+ (r (i+1, j+1) -r (i+1, j) ) 

sr=- ( z (i+1, j+1) -z (i+1, j) ) 

qs=sz*u (ii, j) +sr*v (ii, j) 

qsi=sz*u (i, j)+sr*v(i,  j) 

qsipl=sz*u (i+1, j) +sr*v (i+1, j) 

if  (qsipl . gt . qsi ) qs=0 . 5* (qsi+qsipl ) 

piij=p (ii, j) 

fsl (k2) =rho (ii, j, n) *qs 

fs2  (k2) =rhou (ii, j,n) *qs+piij*sz 

fs3  (k2) =rhov (ii,  j, n) *qs+pii j*sr 

fs4(k2)  = (e(ii,j,n)+piij)*qs 

if  (iflxs.gt.0)  call  fsifs  (i, j, ii) 

if  (nadv. It .nvisc)  go  to  11 

uii j=u (ii, j) 

vii j=v (ii, j ) 

due=u (i+1, j) -u  (i, j) 

dve=v (i+1, j) -v (i, j) 

dee=ei (i+1, j) -ei (i, j) 

dze=0 .25* (z (i+2, j + 1) +z (i+2, j)-z  (i, j + 1) -z (i, j) ) 

dre=0 .25* (r (i+2, j+1) +r (i+2, j) -r (i, j+1) -r (i, j) ) 

dun=0.5* (u (ii, j+1) -u (ii, j— 1 ) ) 

dvn=0 . 5* (v (ii, j+1) -v (ii, j— l ) ) 

den=0 . 5* (ei (ii, j+1) -ei (ii, j-1) ) 

dzn=r (i+1, j+1) -z (i+1, j) 

drn=r (i+1, j+1) -r (i+1, j) 

dzr=1.0/ (dze*drn-dzn*dre) 

duz= (due*drn-dun*dre) *dzr 
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dvz= (dve*drn-dvn*dre) *dzr 
dez= (dee*drn-den*dre) *dzr 
dur=- (due*dzn-dun*dze) *dzr 
dvr=- (dve*dzn-dvn*dze) *dzr 
der=- (dee*dzn-den*dze) *dzr 
temp=abs (0.5*  (ei (i, j) +ei (i  + 1,  j) ) /cv) 
c - use  Sutherland's  form  for  air  in  mks  unit- 
c - rmu=l . 458e-6*sqrt (temp**3) / ( temp+1 1 0 . 4 ) 
c - uf4  viscosity  at  2000  K. 
c rmu=3 . 357e-6*sqrt (temp) / (aa+bb*temp) 

rmu=8 . 67e-5 
c - conductivity, 

c 

rk=gamma*rmu/pr 

c - consider  thermal  radiation  heat  transfer- 
if (nradq.eq. 1) then 

c consider  radiative  conductivity- 

sigma=5 . 67e-8 
tt=ei (i,  j) /cv 

radk= (16 . * sigma* (tt**3))/(3. *rabsco*cv) 
c rabsco-rosseland  absorption  coefficient- 

rk=rk+radk/cv 
endif 
c 
c 

rlmbda=- (2 . 0/3 . 0) *rmu 
stress=-rlmbda* (duz+dvr) 
sigz=stress-2 . 0*rmu*duz 
sigr=stress-2 . 0*rmu*dvr 
tau=-rmu* (dur+dvz) 
qz=rk*dez 
qr=rk*der 
11  continue 
return 
end 



subroutine  fsj(i,j,jj) 



c 

common/a/rho  (100,  100,2), rhou (100,  100,2), rhov(100,  100,2), 
c e (100,  100,2)  , fsl (2)  , fs2 (2)  , 

c u (100,  100) , v (100, 100) , ei (100,  100) , p (100,  100)  , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101, 101) , 

c gamma, gml, ggml, cv, pr,  vou, ptt, ttt , pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl, k2, 

c iadbwl, eiwall, lent , cfl, 

c nadv, nend, nvisc 

common/b/iturb, prt, ep (100,  100) 

common/c/du (100, 100,4), dus (100,100,4) 
common/ e/gl 1 (101 ) , gl2 (101 ) ,gl3 (101),  gl 4 (101), 
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c g21 (101)  ,g22 (101) ,g23 (101) ,g24 (101) , 

c g31 (101) ,g32 (101) , g33 (101) ,g34 (101) , 

c g41 (101)  , g42 (101) , g43 (101) , g44 (101) , 

c f f (101) , implt 

common/ splt/iflxs, sz, sr 
common/ rad/nradq, rabsco 

aa=0 . 80 
bb=-7 . 10e-5 

sz=- (r (i+1, j+1) -r (i, j+1) ) 

sr=+  (z (i+1,  j + 1) -z (i, j + 1) ) 

qs=sz*u (i, j j) + sr*v (i, j j) 

qs j=sz*u (i, j) +sr*v (i, j) 

qs jpl=sz*u (i, j+1) +sr*v (i, j+1) 

if  (qsjpl.gt.qsj)  qs=0 . 5* (qs j+qs jpl) 

pij j=p  (i, j j) 

if  ( j . eq. 1 . or . j . eq. jlml ) qs=0.0 

fsl (k2) =rho (i, j j, n) *qs 

fs2 (k2) =rhou (i, j j,  n) *qs+pijj*sz 

fs3 (k2) =rhov ( i , j j, n) *qs+pi j j*sr 

fs4(k2)  = (e(i,jj,n)+pijj)*qs 

if  (iflxs.gt.0)  call  fs jfs  (i, j, j j) 

if  (nadv. It .nvisc)  go  to  11 

ui j j=0 . 5* (u(i,  j) +u (i, j+1) ) 

vij  j=0.5*  (v(i,  j)  +v(i,  j + 1)  ) 

due=0.5*  (u (i+1,  j j)  -u(i-l,  jj)  ) 

dve=0 . 5* (v (i+1 , j j) — v (i— 1, j j) ) 

dee=0 . 5* (ei (i+1, jj)-ei (i-1, jj) ) 

dze=z (i+1, j+1) -z (i,  j + 1) 

dre=r (i+1, j+1) -r (i,  j + 1) 

dun=u (i, j + 1) -u  (i, j) 

dvn=v (i, j+1) -v (i, j) 

den=ei (i, j+1) -ei (i, j) 

dzn=0 .25* (z (i+1,  j+2)+z (i, j+2)-z  (i  + 1, j)-z(i, j) ) 
drn=0 . 25* (r (i+1, j+2) +r (i, j+2) -r (i+1, j) -r (i, j) ) 
dzr=1.0/ (dze*drn-dzn*dre) 
duz= (due*drn-dun*dre) *dzr 
dvz= (dve*drn-dvn*dre)  *dzr 
dez= (dee*drn-den*dre) *dzr 
dur=- (due*dzn-dun*dze) *dzr 
dvr=- (dve*dzn-dvn*dze) *dzr 
der=- (dee*dzn-den*dze)  *dzr 
temp=abs (0 . 5* (ei (i,  j) +ei (i, j + 1) ) /cv) 
c -Sutherland's  form  for  air  viscosity  in  mks  unit- 
c rmu=l . 458e-06*sqrt (temp* *3) / (temp+110 . 4) 

c - uf 4 viscosity  at  2000  K 
c rmu=3 . 357e-6*sqrt (temp) / (aa+bb*temp) 

rmu=8 . 67e-5 
c 

rk=gamma*rmu/pr 
rmu=rmu+ep (i,  j) 
rk=rk+gamma*ep (i, j) /prt 


c 

c 
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c 

c 


c 


c 

c 

c 

c 


- consider  radiation  heat  transfer- 
if (nradq. eq. 1) then 

- consider  radiative  conductivity- 
sigma=5 . 67e-8 

- rabsco-rosseland  absorption  coefficient- 
tt=ei  (i, j ) /cv 

radk= (16 . * sigma* (tt**3)  ) / (3.  *rabsco*cv) 
rk=rk+radk/cv 
endif 

rlmbda=- (2 . 0/3. 0) *rmu 
stress=-rlmbda* (duz+dvr) 
sigz=stress-2 . 0*rmu*duz 
sigr=stress-2 . 0*rmu*dvr 
tau=-rmu* (dur+dvz) 
qz=rk*dez 
qr=rk*der 

fs2 (k2)=fs2 (k2)+sigz*sz+tau*sr 
fs3 (k2) =fs3  (k2) +tau*sz+sigr*sr 
fs4 (k2) =fs4  (k2)  + (sigz*uijj+tau*vi j j-qz ) *sz  + 
c (tau*ui j j+sigr*vi j j-qr) *sr 

11  continue 
return 
end 


subroutine  fsh(i,j,jj) 


common/a/rho  (100,  100,2) ,rhou(100,  100,2) ,rhov(100,  100,2) , 
c e (100, 100,2) , fsl (2) , fs2  (2) , 

c u (100, 100) , v (100,  100)  ,ei (100,  100) , p (100,  100) , 

c fs3 (2) , fs4 (2) , z (101,  101) , r (101, 101) , vol (101, 101) , 

c gamma, gml, ggml,  cv,pr,  vou,ptt, ttt,pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl, k2, 

c iadbwl, eiwall, lent , cf 1, 

c nadv, nend, nvisc 

common/b/iturb, prt, ep (100, 100) 
common/c/du (100,  100,  4),  dus (100,  100,4) 
common/ e/gll (101),gl2(101),gl3(101),gl4(101), 
c g21 (101) ,g22 (101)  ,g23 (101) ,g24  (101) , 

c g31  (101)  ,g32  (101)  ,g33  (101)  ,g34  (101)  , 

c g41 (101) ,g42 (101)  ,g43 (101) ,g4  4 (101) , 

c f f (101) , implt 

common/spit /if lxs, sz, sr 
common/q/qp 

common/rad/nradq, rabsco 


aa=0 . 80 
bb=-7 . 10e-5 

sz=- (r (i+1, j + 1) -r (i,  j + 1)  ) 
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sr=+ (z (i+1, j+1) -z (i, j+1) ) 

qs=sz  *u ( i , j j) +sr*v (i, j j) 

qs j=sz*u (i, j) +sr*v (i, j) 

qs jpl=sz*u (i, j+l)+sr*v(i,  j + 1) 

if  (qs jpl . gt . qs j)  qs=0 . 5* (qs j+qs jpl) 

Pij j=P (if j j) 

if  ( j . eq. 1 . or . j . eq. jlml)  qs=0 . 0 
fsl (k2) =-rho (i, j j, n) *qs 
fs2 (k2) =-rhou (i, j j,n) *qs 
fs3 (k2) =-rhov (i, j j, n) *qs 
fs4(k2)=-(e(i,jj,n)  +pi  j j)  *qs 
if  (iflxs.gt.O)  call  f s jf s (i, j , j j ) 
if  (nadv. It . nvisc)  go  to  11 
uijj=0.5*  (u  (i, j)+u(i, j + 1) ) 
vi j j=0 .5* (v (i, j) +v (i, j+1) ) 
c 

ry=0 .5* (r (i, j) +r (i, j+1) ) 
due=0 .5* (u  (i  + 1, j j) -u  (i-1, jj) ) 
dve=0 .5* (v (i+1, j j) -v (i-1, jj) ) 
dee=0 . 5* (ei  (i+1, j j) -ei (i-1, j j) ) 
dze=z (i+1, j+1) -z (i, j+1) 
dre=r (i+1, j+1) -r (i, j+1) 
dun=u (i, j+1) -u  (i, j) 
dvn=v (i, j+1) -v (i, j) 
den=ei  (i, j+1) -ei (i, j) 

dzn=0 .25* (z  (i+1, j+2) +z (i, j+2) -z (i+1, j) -z (i, j) ) 
drn=0 . 25* (r (i+1, j+2) +r (i,  j+2) -r (i+1, j) -r (i, j) ) 
dzr=l . 0/ (dze*drn-dzn*dre) 
duz= (due*drn-dun*dre) *dzr 
dvz= (dve*drn-dvn*dre) *dzr 
dez= (dee*drn-den*dre) *dzr 
dur=- (due*dzn-dun*dze) *dzr 
dvr=- (dve*dzn-dvn*dze) *dzr 
der=- (dee*dzn-den*dze) *dzr 
temp=abs  (0.5*(ei(i,j)  +ei  (i,  j + 1)  ) /cv) 
c - uf4  viscosity  at  2000  K 
rmu=8 . 67e-5 

c rmu=3 . 357e-6*sqrt (temp) / (aa+bb*temp) 
c - viscosity  for  air. 

c 

c rmu=l . 458e-6*sqrt (temp* *3)  / (temp+110 . 4 ) 

rk=gamma*rmu/pr 
rmu=rmu+ep ( i , j ) 
rk=rk+gamma*ep (i, j) /prt 
c - consider  radiation  heat  transfer- 

if (nradq. eq. 1) then 

c - consider  radiative  conductivity- 

sigma=5 . 67e-8 

c rabsco-rosseland  absorption  coefficient- 

tt=ei (i, j) /cv 

radk= (16. *sigma* (tt**3)  ) / (3 . *rabsco*cv) 
rk=rk+radk/cv 


endif 
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c 

c 

r lmbda=- ( 2 . 0 / 3 . 0 ) * rmu 
stress=-rlmbda* (duz+dvr) 
sigz=stress-2 . 0*rmu*duz 
sigr=stress-2 . 0*rmu*dvr 
tau=-rmu* (dur+dvz) 
qz=rk*dez 
qr=rk*der 

fs2  (k2)=fs2 (k2) +sigz*sz+tau*sr 
fs3 (k2) =fs3 (k2) +tau*sz+sigr*sr+stress 
c - (4 . /3 . ) *rmu*ui j j/ry 

fs4(k2)=fs4(k2)+ (sigz*ui j j+tau*vi j j-qz) *sz+ 
c (tau*ui j j+sigr*vi j j-qr) *sr 

11  continue 
return 
end 


c 

subroutine  l(dt) 
c 


common/ a/ rho (100,100,2), rhou (100,100,2), rhov (100, 100,2) , 
c e (100, 100, 2) , fsl (2) , fs2 (2) , 

c u (100, 100) , v (100, 100) ,ei(100,100) , p (100, 100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol  (101,101), 

c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl, k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

common /c/du (100,100,4), dus (100,100,4) 
common/e/gll (101),gl2(101),gl3(101),gl4(101), 
c g21 (101) ,g22 (101) ,g23 (101) ,g24 (101) , 

c g31 (101) , g32 (101) , g33 (101)  ,g34 (101) , 

c g41 (101) , g42 (101) ,g43 (101) ,g44 (101) , 

c f f ( 101 ) , implt 

common /spit/ iflxs,sz,sr 
common/tem/t (101, 101)  , ttwall 
common/q/qp 

lcnt=lcnt+l 

iadd=mod (nadv,  2) 

jadd=mod ( ( (iadd+nadv) /2)  , 2) 

do  20  n=l,2 

m=3-n 

nml=n-l 

call  lbc (dt ) 

do  5 j=2, jlml 

call  f si ( 1 , j , 1+iadd) 

do  5 i=2,ilml 

k=kl 
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kl=k2 

k2=k 

call  f si  (i, j , i+iadd) 
dtv=dt/vol  (i,  j) 

du (i, j, 1) =-dtv* (fsl(k2)-fsl(kl)) 
du  (i,  j,  2)  =-dtv*  (fs2(k2)-fs2(kl)) 
du  (i,  j,  3)  =-dtv*  (fs3(k2)-fs3(kl)  ) 
du  (i,  j,  4) =-dtv* (fs4 (k2)-fs4 (kl) ) 

5 continue 

do  10  i=2,ilml 

call  fs j (i, 1, 1+ jadd) 

do  10  j=2, jlml 

k=kl 

kl=k2 

k2=k 

call  fs j (i, j, j+ jadd) 
dtv=dt/vol (i, j) 

du  (i,  j,  1) =du (i, j, 1) -dtv* (fsl (k2)-fsl (kl) ) 
du (i, j, 2) =du (i, j, 2) -dtv* (fs2 (k2)-fs2 (kl) ) 
du  (i,  j,  3)  =du  (i,  j,  3)  -dtv*  (fs3(k2)-fs3(kl)  ) 
du (i, j, 4) =du (i, j, 4) -dtv* (fs4 (k2)-fs4  (kl) ) 

10  continue 
c 

do  11  i = 2, ilml 
call  fsh (i, 1, 1+ jadd) 
do  11  j = 2, jlml 
call  fsh (i, j, j+ jadd) 
dtv=dt/vol (i, j) 

du  (i,  j,  1)  =du  (i,  j,  1)  +dtv*  (fsl(k2)*vol(i,j)) 
du (i,  j, 2) =du (i, j, 2) +dtv* (fs2 (k2) *vol (i, j) ) 
du  (i,  j,  3)  =du (i,  j,  3) +dtv* (fs3 (k2) *vol  (i, j) ) 
du  (i,  j,  4)  =du  (i,  j,  4)  +dtv*  (fs4(k2)*vol(i,j))  +qp*dt 

11  continue 
c 

if  (implt.gt.O)  call  limplt (dt , iadd, jadd) 
rn=n 

rnml=rn-l . 0 
do  12  i=l,il 
do  12  j=2, jlml 

rho (i, j, m) = (rnml*rho (i, j , m) +rho (i, j , n) +du (i, j, 1 ) ) /rn 
rhou (i, j, m)  = (rnml*rhou (i, j,  m) +rhou (i, j , n) +du (i,  j, 2) ) /rn 
rhov (i,  j, m)  = (rnml*rhov (i, j , m) +rhov (i, j , n) +du (i, j, 3) ) /rn 
e (i,  j,m)  = (rnml*e  (i,  j,  m)  +e  (i,  j,  n)  +du  (i,  j,  4) ) /rn 
u (i,  j)  =rhou (i, j, m) /rho  (i, j,m) 
v (i,  j)  =rhov (i, j,m) /rho (i, j, m) 

ei  (i,  j)  =e  (i,  j,m)  /rho  (i,  j,  m)  -0 . 5*  (u  (i,  j)  *u  (i,  j)  +v  (i,  j)  *v  (i,  j)  ) 
p (i, j) =gml*rho (i, j , m) *ei (i , j ) 
c 

t (i, j) =ei  ( i , j ) /cv 
c 

12  continue 
call  be 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


iadd=mod ( (iadd+1 ) , 2) 
j add=mod ( ( j add+ 1 ) , 2 ) 
20  continue 
return 
end 


if  1-D  integral  approximation  is  used,  this  subroutine,  rq, 
must  be  included. 


subroutine  rq 


common/ a/ rho (100,100,2), rhou (100,100,2)  , rhov (100,100,2), 
c e (100,  100, 2) , fsl  (2) , fs2 (2) , 

c u (100, 100) , v (100, 100) ,ei (100, 100) ,p (100, 100) , 

c fs3(2) , f s4 (2 ) , z ( 101 , 1 01 ) , r ( 1 0 1 , 1 01 ) , vol (101, 101)  , 

c gamma, gml, ggml, cv, pr, vou, ptt , ttt, pe, 

c il,jl,ilml,jlml,jlfm,n,m,kl,k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

common/ c/du (100,100,4) , dus (100,100,4) 
common /e/gll (101) ,gl2 (101) ,gl3 (101) ,gl4  (101) , 
c g21 (101) ,g22 (101) ,g23 (101) ,g24 (101) , 

c g31 (101) , g32 (101) , g33 (101) ,g34  (101) , 

c g41 (101) ,g42 (101) ,g43 (101) ,g44 (101) , 

c ff (101) , implt 

common/spit /if lxs, sz, sr 
common/tem/t (101, 101) , ttwall 
common/q/qp 
common/ r/radq (100,100) 
common/rad/nradq,  rabsco 

dimension  rr(100),  plusi (100, 100) , rminusi (100, 100) 

calculate  + direction  radiaitve  intensity 

pi=3. 14159 

sigma=5 . 67e-8 

do  1 i -2 , il-1 

do  2 k =2, j 1-1 

r=0 . 0 

plusin=0 . 0 

do  3 j =2, jl-2 

rr  ( j)  =0 . 5*  (r  (i,  j+2)  -r  (i,  j)  ) 

r=r+rr ( j ) 

if (r . le . (5 . /rabsco) ) then 
ri= ( (t (i,  j) +t (i, j + 1) ) *0.5) **4 

plusi (i, j) =rabsco*sigma*ri* (1 . /pi) *exp ( -rabsco *rr  ( j) ) *rr ( j) 

plusin=plusin+plusi  (i, j) 

endif 


c3 

c 

c2 


continue 

plusi (i, k) =plusin*2. *pi*r (i,k+l) * (z (i  + 1,  k) -z (i, k) ) 
continue 
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c 

c - calculate  - direction  radiative  intensity 
c aaw=rabsco*  (ttt/ttwall) 

c do  4 k= jl, 4 , -1 

c rminusin=0.0 

c r=0.0 

c do  11  j=jl,4,-l 

c rw=r (i, j) -r (i, 2) 

c rr ( j) =0 . 5* (r  (i,  j)  -r  (i, j-2) ) 

c r=r+rr(j) 

c if (r . le . (5 . /rabsco) . and.  rw . gt . (5 . /aaw) ) then 

c ri= ( (t (i, j) +t (i, j-1) ) *0 . 5) **4 

c rminusi (i, j) =rabsco*sigma*ri*  (1 . /pi) * 

c * exp (-rabsco*rr ( j) ) *rr ( j) 
c rminusin=rminusin+rminusi (i, j) 

c endif 

ell  continue 

c rminusi ( i , k) =rminusin 

c if  (r . le . (5 . /rabsco)  . and.  rw . le . (5 . /aaw) then 

c rminusi (i, k) =rminuin+ (1 . /pi) *sigma* (ttwall**4) *exp (-rw*aaw) 

c endif 

c4  continue 

cl  continue 

c - calculate  radiative  energy  deposited  in  every  cell, 
c 

c do  12  i = 2,  il-1 

c do  12  k = 2,  j 1-2 

c radq (i, k) =plusi (i, k) -rminusi (i, k) 

cl2  continue 

c return 

c end 

c 
c 

c 

subroutine  impli (dt, i, j,  iadd) 

c 

c 

c -implicit  procedure  in  i-  direction. 

common/ a/ rho (100,100,2)  , rhou (100,100,2) , rhov (100,100,2), 
c e (100, 100,2) ,fsl (2) , fs2  (2) , 

c u (100, 100) , v (100, 100) ,ei (100, 100) , p (100, 100) , 

c fs3  (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml, ggml , cv,  pr, vou, ptt , ttt , pe, 

c il, jl, ilml , jlml , jlfm, n, m, kl , k2, 

c iadbwl, eiwall,  lent,  cfl, 

c nadv, nend,  nvisc 

common /b/ iturb, prt, ep (100, 100) 
common/ c/du (100,100,4), dus (100, 100, 4) 
common/d/all, al2,al3,al4,a21,a22,a23,a24,ffi(2), 
c a31,a32,a33,a34,a41,a42,a43,a44,ni 

common/e/ gll (101) ,gl2 (101),  gl3 (101) ,gl 4(101), 
c g21 (101) ,g22 (101) ,g23 (101) , g24  (101) , 
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c 

c 

c 


g31 (101) ,g32 (101) ,g33 (101) ,g34 (101) , 
g41 (101) , g42 (101) , g43 (101) ,g44 (101) , 
ff (101) , implt 


common/splt/if lxs, sz, sr 


c 

c 


beta=gamma-l . 0 
sign=l . 0 
do  20  nn=l,2 
ii=i+2-nn 

si=sqrt ( (z (ii, j + l)-z (ii, j) ) **2+ (r  (ii,  j + l)-r(ii, j) ) **2) 
szp=  (r (ii, j+1) -r (ii, j) ) /si 
srp=- (z (ii, j+1) -z (ii, j) ) /si 
ii=ii-l+iadd 

c if  (ii.eq.il)  ii=ii-l 

rr=p (ii, j) / (gml*ei (ii, j) ) 
uu=u (ii, j) 
vv=v (ii, j) 

cc=sqrt (ggml*abs (ei (ii, j) ) ) 
up=  szp*uu+srp*vv 
vp=-srp*uu+szp*vv 
alp=0 . 5* (uu*uu+vv*vv) 
rccsq=l .0/ (cc*cc) 
dtvol=dt/vol (i,  j) 
if lag=l 

if  (ni . eq. 2 . and. nn . eq. 1)  iflag=2 
6 continue 

il=i+2-nn+l-if lag 
i2=i 1-3+2 *iflag 
aa=l . 0 

if  (iflxs.eq.2)  aa=l . 5 
dd=-aa*sign* (3-2* (iflag) ) *si*dtvol 
dd0=0 . 0 

if  (abs  (p(il,  j) -p(i2,  j)  ) .gt. 
c 0 . 5*aminl (abs (p (il, j) ) , abs (p (i2, j) ) ) ) then 
rr=rho (il , j , n) 
uu=u  (il, j) 
vv=v (il, j) 
up=  szp*uu+srp*vv 
vp=-srp*uu+szp*vv 
cc=sqrt (ggml*abs (ei (il, j) ) ) 
alp=0 . 5* (up*up+vp*vp) 

dd0=0 . 5*abs (dd) *amaxl (cc-abs (up) ,0.0) /gamma 
endif 

ddl=up*dd 
dd2= (up+cc) *dd 
dd4= (up-cc) *dd 
il=i+2-nn 

coef=l . 0 

if  (il . eq. 2 . or . il . eq. il)  coef=0.0 
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epl=coef *dd*  (szp* (u(il, j) -u (il-1, j) ) +srp* (v (il, j) 
ep2=coef *dd* (szp* (u (il, j) -u (il-1, j) ) +srp* (v(il, j 
c + (sqrt  (ggml*abs  (ei  (il,  j)  ) ) -sqrt  (ggml*abs  (ei  (il- 
ep4=coef*dd* (sxp*(u(il,j)-u(il-l,j))+syp*(v(il,j 
c - (sqrt  (ggml*abs  (ei  (il,  j)  ) ) -sqrt  (ggml*abs  (ei  (il- 

ddl=0 . 5* (ddl+sqrt (ddl *ddl+epl *epl ) ) +ddO 
dd2=0 .5* (dd2+sqrt (dd2*dd2+ep2*ep2) ) +ddO 
dd4=0 .5* (dd4+sqrt (dd4*dd4+ep4*ep4 ) ) +ddO 

il=i+2-nn+l-if lag 
tll=  ddl* (1 . 0-alp*beta*rccsq) 
tl2=  ddl*uu*beta*rccsq 
tl3=  ddl*vv*beta*rccsq 
t 14=-ddl*beta*rccsq 
t21=  dd2* (alp*beta-up*cc) 
t22=  dd2* ( sxp*cc-uu*beta) 
t23=  dd2* ( syp*cc-vv*beta) 
t24=  dd2*beta 
t31=-ddl*vp/rr 
t32=-ddl*syp/rr 
t33=  ddl*sxp/rr 
t 34=  0.0 

t41=  dd4* (alp*beta+up*cc) 
t42=  dd4*  (-szp*cc-uu*beta) 
t43=  dd4* (-srp*cc-vv*beta) 
t44=  dd4*beta 

if  (iflag.eq.2)  go  to  8 
bpll=  1.0 
bpl2=  0.5*rccsq 
bpl3=  0.0 
bpl4=  0.5*rccsq 
bp21=  uu 

bp22=  0 . 5*  (uu+szp*cc) *rccsq 
bp23=-srp*rr 

bp24=  0.5*  (uu-szp*cc) *rccsq 
bp31=  vv 

bp32=  0 . 5* (vv+srp*cc) *rccsq 
bp33=  szp*rr 

bp34=  0 . 5* (vv-srp*cc) *rccsq 
bp41=  alp 

bp42=  0 . 5*  ( (alp+up*cc) *rccsq+l . 0/beta) 
bp43=  rr*vp 

bp4  4=  0 . 5*  ( ( alp-up *cc) *rccsq+l . 0/beta) 
bll=- (bpll*tll+bpl2*t21+bpl3*t31+bpl4*t41) 
bl2=- (bpl 1 *t 12+bpl2*t22+bpl3*t32+bpl4*t 42 ) 
bl3=- (bpl 1 *t 13+bpl2*t23+bpl3*t33+bpl4  *t43) 
bl4=- (bpll*t 14+bpl2*t24+bpl3*t34+bpl4*t44 ) 
b21=- (bp21*tll+bp22*t21+bp23*t31+bp24*t41) 
b22=- (bp21 *t 12+bp22*t22+bp23*t 32+bp24  *t 42 ) 
b23=- (bp21*tl3+bp22*t23+bp23*t33+bp24*t43) 


■v (il-1, j) ) ) 
-v (il-1, j)  ) 
1,  j) ) ) ) ) 

-v (il-1, j) ) 
1, j) ) ) ) ) 
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b24=- (bp21*tl4+bp22*t24+bp23*t34+bp24*t44 ) 
b31=- (bp31*tll+bp32*t21+bp33*t31+bp34*t41) 
b32=- (bp31*t 12+bp32*t22+bp33*t32+bp34  *t42 ) 
b33=- (bp31*tl3+bp32*t23+bp33*t33+bp34*t43) 
b34=- (bp31*t 14+bp32*t24+bp33*t34+bp34*t44 ) 
b4 1=- (bp41*tll+bp42*t21+bp43*t31+bp44*t41) 
b42=- (bp41*t 12+bp42*t22+bp4  3*t32+bp4  4 *t 42 ) 
b43=- (bp41*t 13+bp42*t23+bp43*t33+bp44  *t43) 
b44=- (bp4 1 *t 14+bp42*t24+bp4  3*t 34+bp44  *t 4 4 ) 
if lag=2 
ab=l . 0 
ac=0 . 0 

if  (nn.eq.l)  go  to  12 
go  to  14 
8 continue 
cp 11=  1.0 
cpl2=  0.5*rccsq 
cpl3=  0.0 
cpl4=  0.5*rccsq 
cp21=  uu 

cp22=  0 . 5* (uu+szp*cc) *rccsq 
cp23=-syp*rr 

cp24=  0 . 5* (uu-szp*cc) *rccsq 
cp31=  vv 

cp32=  0 . 5* (vv+srp*cc) *rccsq 
cp33=  szp*rr 

cp34=  0 . 5* (vv-srp*cc) *rccsq 
cp41=  alp 

cp42=  0 . 5* ( (alp+up*cc) *rccsq+l . 0/beta) 
cp43=  rr*vp 

cp44=  0 . 5* ( ( alp-up *cc) *rccsq+l . 0/beta) 
cll=- (cpll*tll+cpl2*t21+cpl3*t31+cpl4*t41) 
cl2=- (cpll*tl2+cpl2*t22+cpl3*t32+cpl4*t42) 
cl3=- (cpll*tl3+cpl2*t23+cpl3*t33+cpl4*t43) 
Cl4=- (cpll*t 14+cpl2*t24+cpl3*t34+cpl4*t44 ) 
c21=- (cp21 *t Il+cp22*t21+cp23*t 31+cp24  *t4 1 ) 
c22=- (cp21*t 12+cp22*t22+cp23*t32+cp24*t 42) 
c23=- (cp21*t 13+cp22*t23+cp23*t33+cp24*t 43) 
c24=- (cp21*t 14  + cp22*t24+cp23*t34+cp24  *t 44 ) 
c31=- (cp31*t Il+cp32*t21+cp33*t31+cp34*t41 ) 
c32=- (cp31*tl2+cp32*t22+cp33*t32+cp34*t42) 
c33=- (cp31*t 13+cp32*t23+cp33*t 33+cp34*t 43) 
c34=- (cp31*t 14+cp32*t24+cp33*t34+cp34*t44) 
c4 1=- (cp41*tll+cp42*t21+cp43*t31+cp44*t41) 
c42=- (cp41*tl2+cp42*t22+cp43*t32+cp44*t42) 
c43=- (cp41*t 13+cp42*t23+cp43*t33  + cp4  4 *t 43) 
044=-(cp41*tl4+cp42*t24+cp43*t34+cp44*t44) 
ab=0 . 0 
ac=l . 0 

if  (nn.eq.l)  go  to  14 
if  (ni.eq.l)  go  to  20 
il=i-l 
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dul=du (il, j, 1) 
du2=du (il, j, 2) 
du3=du (il , j , 3) 
du4=du (il, j, 4) 

du (i, j , 1 ) =du (i, j, l)-(cll*dul+cl2*du2+cl3*du3+cl4*du4) 
du (i, j, 2) =du (i, j, 2) - (c21*dul+c22*du2+c23*du3+c24*du4) 
du (i, j, 3) =du (i, j, 3) - (c31*dul+c32*du2+c33*du3+c34*du4) 
du (i, j, 4) =du (i, j, 4) - (c41*dul+c42*du2+c43*du3+c44*du4) 
go  to  20 
12  continue 
il=i+l 

dul=du (il, j, 1) 
du2=du (il, j, 2) 
du3=du (il, j, 3) 
du4=du (il, j, 4) 

bdul=bll*dul+bl2*du2+bl3*du3+bl4*du4 
bdu2=b21*dul+b22*du2+b23*du3+b24  *du4 
bdu3=b31 *dul+b32*du2+b33*du3+b34  *du4 
bdu4=b41*dul+b42*du2+b43*du3+b44*du4 
du ( i , j , 1 ) =du ( i , j , 1 ) -bdul 
du ( i , j , 2 ) =du ( i , j , 2 ) -bdu2 
du ( i , j , 3 ) =du ( i , j , 3 ) -bdu3 
du ( i , j , 4 ) =du ( i , j , 4 ) -bdu4 
dus ( i , j , 1 ) =dus ( i , j , 1 ) -bdul 
dus ( i , j , 2 ) =dus ( i , j , 2 ) -bdu2 
dus ( i , j , 3 ) =dus ( i , j , 3 ) -bdu3 
dus (i, j, 4) =dus (i, j, 4) -bdu4 
go  to  6 
14  continue 

all=all-ab*bll-ac*cll 
al2=al2-ab*bl2-ac*cl2 
al3=al3-ab*bl3-ac*cl3 
al4=al4-ab*bl4-ac*cl4 
a21=a21-ab*b21-ac*c21 
a22=a22-ab*b22-ac*c22 
a23=a23-ab*b23-ac*c23 
a24=a24-ab*b24-ac*c24 
a31=a31-ab*b31-ac*c31 
a32=a32-ab*b32-ac*c32 
a33=a33-ab*b33-ac*c33 
a34=a34-ab*b34-ac*c34 
a41=a41-ab*b41-ac*c41 
a42=a42-ab*b42-ac*c42 
a43=a43-ab*b43-ac*c43 
a44=a44-ab*b44-ac*c44 
if  (nn.eq.2)  go  to  6 
20  continue 
return 
end 
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subroutine  limplt (dt , iadd, jadd) 

c 

c 

c****  implicit  procedure  in  j-  direction. 

common/ a/ rho (100,100,2) , rhou (100,100,2),  rhov (100,100,2), 
c e (100,  100, 2) , fsl  (2) , fs2  (2) , 

c u(100, 100) , v (100, 100) ,ei (100, 100) , p (100, 100) , 

c fs3 (2) , fs4  (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm, n, m,  kl,  k2, 

c iadbwl, eiwall, lent,  cfl, 

c nadv, nend, nvisc 

common/b/iturb, prt , ep (100,  100) 
common/ c/du (100,  100,  4)  , dus (100,  100,  4) 
common/d/ all,al2,al3,al4,a21,a22,a23,a24,  ffi (2) , 
c a31,a32,a33,a34,a41,a42,a43,a44,ni 

common/e/gll (101) ,gl2 (101),gl3(101),gl4(101), 
c g21 (101) , g22 (101) ,g23 (101) , g24 (101) , 

c g31 (101) , g32 (101) ,g33 (101) ,g34 (101) , 

c g41 (101) ,g42 (101) ,g43 (101) ,g44 (101) , 

c f f (101 ) , implt 

common/spit /if lxs, sz,  sr 
dimension  dui(100,4) 

data  all,  al2,  al3,  al4,  a21,  a22,  a23,  a24 , a31 , a32 , a33,  a34, 
c a41,a42,a43,a44,bll,bl2,bl3,bl4,b21,b22,b23,b24, 

c b31,b32,b33,b34,b41,b42,b43,b44,cll,cl2,cl3,cl4, 

c c21,c22,c23, c24, c31, c32, c33, c34,  c41,  c42,  c43,  c44/4 8*0.0/ 
c 

beta=  gamma-1.0 
sign=l . 0 
aa=0 . 80 
bb=-7 . 10e-5 
bc=l. 0/3.0 
do  1 i=l , il 
do  1 j = 1 , j 1 
do  1 1=1,4 

dus  (i,  j,  1)  =du  (i,  j,  1) 

1 continue 

c do  20  ni=l,3 

do  20  ni=l,2 
do  20  iz=2, ilml 
i=ilml+2-iz 
if  (ni.eq.2)  i=iz 
if  (ni.eq.l)  go  to  3 
do  2 j=l, jl 
do  2 1=1,  4 
dui  ( j,  1)  =du  (i,  j,  1) 
du  (i,  j,  1)  =dus  (i,  j,  1) 

2 continue 

3 continue 
ttu=l . 0 


vol (i, jl+1 ) =vol (i, jl) 
do  10  jz=2, jl 
j=jl+2-jz 

s j=sqrt  ((z(i+l,j)-z(i,j))**2+(r(i+l,j)-r(i,j))**2) 

srp=  (z (i+1, j) -z (i, j) ) /s j 

szp=-  (r(i  + l,j)-r(i,j))/sj 

vr=vol (i, j + 1) /vol (i, j) 

temp=abs (0.5*  (ei  (i, j) +ei (i, j— 1 ) ) /cv) 

- use  Sutherland's  form  for  air  viscosity  - 
rmu=l . 458e-06*sqrt (temp* *3) / (temp+110 . 4) 

- uf4  property 

rmu=3 . 367e-6*sqrt (temp) / (aa+bb*temp) 

UF 4 viscosity  at  2000  K 
rmu=8 . 67e-5 


if  (nadv . It . nvisc)  rmu=0.0 
rk=gamma*rmu/pr 
rmu=rmu+ep (i, j-1) 
rk=rk+gamma*ep (i, j-1) /prt 
rlmbda=- (2 . 0/3 . 0) *rmu 
j j=j~l+jadd 

rr=p (i, j j) / (gml*ei  (i, j j) ) 
cc=sqrt (ggml*abs (ei (i, j j) ) ) 
rccsq=l . 0/ (cc*cc) 

upbar=0 .5* (srp* (u  (i,  j)  +u  (i,  j-1)  ) -szp*  (v  (i,  j)  +v  (i,  j-1)  ) ) 
vpbar=0 .5*(szp*(u(i,j)+u(i,j-l))+srp*(v(i,j)+v(i,j-l))) 
dtvol=dt/vol (i, j) 
if lag=l 

if  ( j . eq. jl)  go  to  6 
dul=du (i, j+1, 1) 
du2=du (i, j+1, 2) 
du3=du (i, j+1 , 3) 
du4=du (i, j+1, 4) 

du  (i, j, 1) =du  (i, j, 1) -vr* (bll*dul+bl2*du2+bl3*du3+bl4*du4) 
du (i, j, 2) =du (i, j, 2) -vr* (b21*dul+b22*du2+b23*du3+b24*du4) 
du (i, j, 3) =du (i, j, 3) -vr* (b31*dul+b32*du2+b33*du3+b34*du4) 
du (i, j, 4) =du (i, j, 4) -vr* (b41*dul+b42*du2+b43*du3+b44*du4) 
all=bll*gll ( j+1) +bl2*g21 (j+1) +bl3*g31 ( j+1) +bl4*g41 (j+1) 
al2=bl 1 *gl2 ( j+1 ) +bl2*g22 (j+1) +bl3*g32 ( j+1) +bl4*g42 (j+1) 
al3=bll*gl3  ( j + 1) +bl2*g23 ( j + 1) +bl3*g33 ( j + 1) +bl4*g43 ( j + 1) 
al4=bll*gl4  ( j + 1 ) +bl2*g24 (j  + 1) +bl3*g34 ( j + 1 ) +bl4 *g4  4 (j  + 1) 
a21=b21*gll  (j  + 1) +b22*g21 (j  + 1) +b23*g31 (j  + 1) +b24*g41 (j  + 1) 
a22=b2 1 *gl2  (j  + 1) +b22*g22 ( j + 1 ) +b23*g32 (j  + 1) +b24*g42 (j  + 1) 
a23=b21*gl3 (j+1) +b22*g23 (j+1) +b23*g33 (j+1) +b24*g43 (j+1) 
a24=b21*gl4  ( j + 1) +b22*g24 (j  + 1) +b23*g34  ( j + 1) +b24*g44 (j  + 1) 
a31=b31*gll ( j+1) +b32*g21 ( j+1) +b33*g31 ( j+1) +b34*g41 (j+1) 
a32=b31*gl2  ( j + 1 ) +b32*g22 ( j + 1) +b33*g32  ( j + 1) +b34*g42 ( j + 1) 
a33=b31*gl3  ( j + 1 ) +b32*g23 (j  + 1) +b33*g33 ( j + 1) +b34*g4  3 (j  + 1) 
a34=b31 *gl 4 ( j + 1) +b32*g24 ( j + 1 ) +b33*g34  ( j + 1 ) +b34 *g4  4 (j  + 1) 
a41=b41*gll  ( j + 1 ) +b42*g21 ( j + 1 ) +b4  3*g31 ( j + 1 ) +b4  4 *g4 1 (j  + 1) 
a42=b4 1 *gl2 ( j + 1 ) +b42*g22 ( j + 1) +b4  3*g32  (j  + 1) +b44*g42 (j  + 1) 
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a43=b41*gl3 ( j + 1 ) +b42*g23 ( j + 1 ) +b43*g33 ( j + 1 ) +b44*g43 ( j+1) 
a44=b41*gl4 ( j+1 ) +b42*g24 (j  + 1) +b43*g34 (j  + 1) +b44*g44  (j  + 1) 
6 continue 
uu=u  (i, j j) 
w=v(i,  j j) 
up=  srp*uu-szp*vv 
vp=  szp*uu+srp*vv 
alp=0.5* (up*up+vp*vp) 
j 1= j+l-if lag 

rr j=p ( i , j 1 ) / (gml *ei  (i, jl) ) 
uu j=u (i, jl) 
vv  j = v ( i , jl) 
up j=srp*uu j-szp*vv j 
vp j=szp*uu j+srp*vv j 
cc j=sqrt (ggml*abs (ei  (i, jl) ) ) 
alp j=0 . 5* (uu j*uu j+vv j *vv j ) 
aa=l . 0 

if  (iflxs.eq.2)  aa=1.5 
dd=-aa*sign* (3-2* (iflag) ) *sj*dtvol 
dd0=0 . 0 

if  (abs  (p(i,  j)-p(i,  j-i)  ) .gt. 
c 0.5*aminl(abs(p(i,  j)  ) ,abs(p(i,  j-1)  ) ) ) then 
rr=rho (i, jl, n) 
uu=u  (i, jl) 
vv=v (i, jl) 
up=  srp*uu-szp*vv 
vp=  szp*uu+srp*vv 
cc=sqrt (ggml*abs (ei  (i, jl) ) ) 
alp=0 . 5* (up*up+vp*vp) 

dd0=0 . 5*abs (dd) *amaxl (cc-abs (vp) ,0.0) /gamma 

endif 

ttb=l . 0 

ttc=l . 0 

tt=0 . 0 

if  (j.eq.2)  then 

if  (iflag. eq. 2)  go  to  9 

dd=0 . 0 

ttc=0 . 0 

ttu=l . 0 

tte=l-2*iadbwl 
tt=l . 0 
endif 

if  (j.eq.jl)  then 

dd=0 . 0 

ttb=0 . 0 

ttu=-l . 0 

tte=-l . 0 

tt=l . 0 

endif 

ddl=vp*dd 

dd3= (vp+cc) *dd 

dd4= (vp-cc) *dd 
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coef=l . 0 

epl=coef*dd*  (szp*  (u  (i,  j) -u (i,  j-1)  ) +srp*  (v  (i,  j) -v  (i,  j-1)  ) ) 
ep3=coef *dd*  (szp*  (u(i,  j) -u  (i,  j-1)  ) +srp*  (v(i,  j) -v  (i,  j-1)  ) 
c + (sqrt (ggml*abs (ei (i,  j)  ) ) -sqrt (ggml*abs (ei(i,j-l))))) 

ep4=coef*dd*  (szp*  (u(i,  j)  -u  (i,  j-1)  ) +srp*  (v  (i,  j) -v  (i,  j-1)  ) 
c - (sqrt  (ggml*abs  (ei  (i,  j)  ) ) - sqrt  (ggml*abs  (ei(i,j-l))))) 

ddl=0 .5* (ddl+sqrt (ddl *ddl+epl *epl ) ) +ddO 
dd3=0 .5* (dd3+sqrt (dd3*dd3+ep3*ep3) ) +ddO 
dd4=0 . 5* (dd4+sqrt (dd4*dd4+ep4*ep4) ) +ddO 
if  (j.le.jlfm)  then 

wtl=( (z(i,j)-z(i,2))**2+(r(i,j)-r(i,2))**2)/ 
c ( (z (i,  jlfm) -z (i, 2) ) **2+ (r (i, jlfm) -r (i, 2) ) **2) 
wt2=l . O-wt 1 
uu=wtl*uu+wt2*uu j 
vv=wt 1 * vv+wt2  * vv j 
alp=0 . 5* (uu*uu+vv*vv) 
endif 


tll=  ddl*  (1 . 0-alp*beta*rccsq) 
tl2=  ddl*uu*beta*rccsq 
tl3=  ddl *vv*beta*rccsq 
t 14=-ddl *beta*rccsq 
t21=-ddl*up/rr 
t22=  ddl*srp/rr 
t23=-ddl*szp/rr 
t24=  0.0 

t31=  dd3* (alp*beta-vp*cc) 
t32=  dd3* ( szp*cc-uu*beta) 
t33=  dd3* ( srp*cc-vv*beta) 
t 34=  dd3*beta 
t41=  dd4* (alp*beta+vp*cc) 
t42=  dd4* (-szp*cc-uu*beta) 
t43=  dd4* (-srp*cc-vv*beta) 

1 4 4=  dd4*beta 

if  (j.le.jlfm)  then 
vp=wt l*vp+wt2*vp j 
uu=  srp*up+szp*vp 
vv=-szp*up+srp*vp 
alp=0 . 5* (uu*uu+vv*vv) 
endif 

if  (iflag.eq.2)  go  to  8 
bpll=  1.0 
bpl2=  0.0 
bpl3=  0.5*rccsq 
bpl4=  0.5*rccsq 
bp21=  uu 
bp22=  srp*rr 

bp23=  (uu+szp*cc) *0 . 5*rccsq 
bp24=  (uu-szp*cc) *0 . 5*rccsq 
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bp31=  vv 
bp32=-szp*rr 

bp33=  (vv+srp*cc) *0 . 5*rccsq 
bp34=  (vv-srp*cc) *0 . 5*rccsq 
bp41=  alp 
bp42=  rr*up 

bp43=  0 . 5* ( (alp+vp*cc) *rccsq+l . 0/beta) 

bp44=  0 . 5* ( (alp-vp*cc) *rccsq+l . 0/beta) 

vdd=s j*dtvol*s j/ ( . 5*rr j* (vol (i, j) +vol (i,  j— 1 ) ) ) 

rnul=ttb* (1 . 0+tt*ttu) *rmu 

rnu2=ttb* (1 . 0+tt) * (rlmbda+2 . 0*rmu) 

rnu3=ttb* (1 . 0+tt*tte) *rk 

bll=- (bpll*tll+bpl2*t21+bpl3*t31+bpl4*t41) 
bl2=- (bpll*t 12+bpl2*t22+bpl3*t32+bpl4*t42 ) 
bl3=- (bpll*t 13+bpl2*t23+bpl3*t33+bpl4*t43) 
bl4=- (bpll*t 14+bpl2*t24+bpl3*t34+bpl4*t44 ) 
b21=- (bp21*tll+bp22*t21+bp23*t31+bp24*t41) 
c +vdd* ( srp*up j *rnul+szp*vp j *rnu2 ) 

b22=- (bp21*tl2+bp22*t22+bp23*t32+bp24*t42) 
c -vdd* (srp*srp*rnul+szp*szp*rnu2) 

b23=- (bp21*tl3+bp22*t23+bp23*t33+bp24*t43) 
c +vdd*szp*srp* (rnul-rnu2) 

b24=- (bp21*t 14+bp22*t24+bp23*t34+bp24  *t44 ) 
b31=- (bp31*t Il+bp32*t21+bp33*t31+bp34*t41 ) 
c +vdd* (-szp*up j *rnul+srp*vp j*rnu2) 

b32=- (bp31*t 12+bp32*t22+bp33*t32+bp34*t42) 
c +vdd*szp*srp* (rnul-rnu2) 

b33=- (bp31*t 13+bp32*t23+bp33*t33+bp34*t43) 
c -vdd* (szp*szp*rnul+srp*srp*rnu2) 

b34=- (bp31 *tl4+bp32*t24+bp33*t34+bp34*t44 ) 
b41=- (bp41*tll+bp42*t21+bp43*t31+bp44*t41) 
c +vdd* (upbar*up j*rnul+vpbar*vp j*rnu2 

c - (alpj-ei (i, j) ) *rnu3) 

b42=- (bp4 1 *t 12+bp42  *t22+bp4  3*t 32+bp4  4 *t 42 ) 
c -vdd* ( upbar*srp*rnul+vpbar*szp*rnu2-uu j*rnu3) 

b43=- (bp41*tl3+bp42*t23+bp43*t33+bp44*t43) 
c -vdd* (-upbar*szp*rnul+vpbar*srp*rnu2-vvj*rnu3) 

b44=- (bp41*tl4+bp42*t24+bp43*t34+bp44*t44 ) 
c -vdd*rnu3 

if  (j.eq.2)  then 

dl=aa* (1 . 0-vpj/cc j) *s j*dt/vol ( i , j ) 

bpl= (alp j *beta+vp j *cc j ) *dl 

bp2= (-uu j *beta-szp*cc j ) *dl 

bp3= (-vvj*beta-srp*cc j) *dl 

bp4=beta*dl 

b21=b21+szp*bpl 

b22=b22+szp*bp2 

b23=b23+szp*bp3 

b24=b24+szp*bp4 

b31=b31+srp*bpl 

b32=b32+srp*bp2 

b33=b33+srp*bp3 


b34=b34+srp*bp4 

endif 

all=l . -bll-vr* (cll+all) 
al2=  -bl2-vr* (cl2+al2) 
al3=  -bl3-vr*  (cl3+al3) 
al4=  -bl4-vr* (cl4+al4 ) 
a21=  -b21-vr* (c21+a21 ) 

a22=l . -b22-vr* (c22+a22) 
a23=  -b23-vr* (c23+a23) 

a24=  -b24-vr* (c24+a24 ) 

a31=  -b31-vr* (c31+a31 ) 

a32=  -b32-vr* (c32+a32) 

a33=l .-b33-vr*  (c33+a33) 
a34=  -b34-vr* (c34+a34 ) 

a41=  -b4 1-vr * (c4 l+a4 1 ) 

a42=  -b42-vr* (c42+a42) 

a43=  -b43-vr* (c43+a43) 

a44=l.-b44-vr* (c44+a44) 
if lag=2 
go  to  6 
8 continue 
cpll=  1.0 
cpl2=  0.0 
cpl3=  0.5*rccsq 
Cpl4=  0.5*rccsq 
cp21=  uu 
cp22=  srp*rr 

cp23=  (uu+szp*cc) *0 . 5*rccsq 
cp24=  (uu-szp*cc) *0 . 5*rccsq 
cp31=  vv 
cp32=-szp*rr 

cp33=  (vv+srp*cc) *0 . 5*rccsq 
cp34=  (vv-srp*cc) *0 . 5*rccsq 
cp4 1=  alp 
cp42=  rr*up 

cp43=  0 . 5* ( (alp+vp*cc) *rccsq+l . 0/beta) 
cp44=  0 . 5* ( (alp-vp*cc) *rccsq+l . 0/beta) 
vdd=vdd* (p (i, j) / (gml*ei (i,  j)  ) ) /rr j 
rnul=ttc*  (1 . 0+tt *ttu) *rmu 
rnu2=ttc* (1 . 0+tt) * (rlmbda+2 . 0*rmu) 
rnu3=ttc*  (1.0+tt*tte) *rk 

cll=- (cpll*tll+cpl2*t21+cpl3*t31+cpl4*t41) 
cl2=- (cpll*tl2+cpl2*t22+cpl3*t32+cpl4*t42) 
cl3=- (Cpll*tl3+cpl2*t23+cpl3*t33+cpl4*t43) 
cl4=- (cpll *t 14+cpl2*t24+cpl3*t34+cpl4*t44 ) 
c21=- (cp21 *t Il+cp22*t21+cp23*t31+cp24*t41) 
c +vdd* ( srp*up j*rnul+szp*vp j*rnu2) 

c22=- (cp21*t 12+cp22*t22+cp23*t32+cp24*t42) 
c -vdd* (srp*srp*rnul+szp*szp*rnu2) 

c23=- (Cp21*tl3+cp22*t23+cp23*t33+cp24*t43) 
c +vdd*szp*srp* (rnul-rnu2) 

c24=- (cp21*t 14  + cp22*t24+cp23*t34  + cp24  *t44 ) 
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c31=- (Cp31*tll+cp32*t21+cp33*t31+cp34*t41) 
c +vdd* (-szp*up j *rnul+srp*vp j*rnu2) 

c32=- (cp31*t 12+cp32*t22+cp33*t32+cp34*t42) 
c +vdd*szp*srp* (rnul-rnu2) 

c33=- (cp31*tl3+cp32*t23+cp33*t33+cp34*t43) 
c -vdd* (szp*szp*rnul+srp*srp*rnu2) 

c34=- (cp31*tl4+cp32*t24+cp33*t34+cp34*t44) 
c41=- (cp41*tll+cp42*t21+cp43*t31+cp44*t41) 
c +vdd* (upbar*up j*rnul+vpbar*vp j*rnu2 

- (alpj-ei (i, j-1) ) *rnu3) 
c42=- (cp41*tl2+cp42*t22+cp43*t32+cp44*t42) 
c -vdd* ( upbar*srp*rnul+vpbar*szp*rnu2-uu j *rnu3) 

c43=- (cp41*tl3+cp42*t23+cp43*t33+cp44*t43) 
c -vdd* (-upbar*szp*rnul+vpbar*srp*rnu2-vv j*rnu3) 

c44=- (cp41*t 14+cp42*t24+cp43*t34+cp44*t44 ) 
c -vdd*rnu3 

if  (j.eq.jl)  then 

dl=-aa* (1 . O+vpj/ ccj) *sj*dt/vol (i,  j) 

cpl= (alpj*beta) *dl 

cp2= (-uu j*beta) *dl 

cp3= (-vvj*beta) *dl 

cp4=beta*dl 

c21=c21+szp*cpl 

c22=c22+szp*cp2 

c23=c23+szp*cp3 

c24=c24+szp*cp4 

c31=c31+srp*cpl 

c32=c32+srp*cp2 

c33=c33+srp*cp3 

c34=c34+srp*cp4 

go  to  10 

endif 

9 continue 

if  ( j .eq. jl)  go  to  10 

call  impli (dt, i, j, iadd) 

rll=all 

r21=a21 

r31=a31 

r41=a41 

sl2=al2/rll 

sl3=al3/rll 

sl4=al4/rll 

r22=a22-r21*sl2 

r32=a32-r31*sl2 

r42=a42-r41*sl2 

s23= (a23-r21*sl3) /r22 

s24=(a24-r21*sl4) / r22 

r33=a33-r31*sl3-r32*s23 

r43=a43-r41*sl3-r42*s23 

s34=(a34-r31*sl4-r32*s24) /r33 

r44=a44-r41*sl4-r42*s24-r43*s34 

yyl=du (i, j, 1) /rll 


yy2= (du  (i,  j,  2) -r21*yyl) /r22 
yy3= (du(i, j,3)-r31*yyl-r32*yy2) /r33 
yy4  = (du  (i,  j,  4) -r41*yyl-r42*yy2-r43*yy3) /r44 
yy3=yy3-s34*yy4 
yy2=yy2-s23*yy3-s24*yy4 
yyl=yyl-sl2*yy2-sl3*yy3-sl4*yy4 
du (i, j, 1) =yyl 
du (i,  j, 2) =yy2 
du (i, j, 3) =yy3 
du (i, j, 4) =yy4 
if  (j.eq.2)  go  to  10 
yyl=cll/rll 
yy2= (c21-r21*yyl) /r22 
yy3= (c31-r31*yyl-r32*yy2) /r33 
g41  ( j)  = (c41-r41*yyl-r42*yy2-r43*yy3) /r44 
g31  ( j ) =yy3-s34  *g4 1 ( j) 
g21 ( j) =yy2-s23*g31 ( j) -s24*g41 ( j) 
gll ( j) =yyl-sl2*g21 (j)-sl3*g31 (j)-sl4*g41 (j) 
yyl=cl2/rll 
yy2= (c22-r21 *yyl ) /r22 
yy3= (c32-r31*yyl-r32*yy2) /r33 
g42 ( j) = (c42-r41*yyl-r42*yy2-r43*yy3) /r44 
g32  ( j) =yy3-s34*g42 ( j) 
g22 ( j ) =yy2-s23*g32 ( j) -s24*g42 ( j) 
gl2 ( j ) =yyl-sl2*g22 (j)-sl3*g32 ( j ) -sl4 *g42 ( j ) 
yyl=cl3/rll 
yy2= (c23-r21*yyl) /r22 
yy3= (c33-r31*yyl-r32*yy2) /r33 
g43 ( j) = (c43-r41*yyl-r42*yy2-r43*yy3) /r44 
g33 ( j) =yy3-s34*g43 ( j) 
g23 ( j) =yy2-s23*g33 ( j) -s24*g43 ( j) 
gl3 ( j)=yyl-sl2*g23 ( j)-sl3*g33 ( j)  -sl4*g43 ( j) 
yyl=cl4/rll 
yy2= (c24-r21*yyl) /r22 
yy3= (c34-r31*yyl-r32*yy2) /r33 
g4  4 ( j)  = (c4  4-r41*yyl-r42*yy2-r4  3*yy3) /r4  4 
g34 ( j) =yy3-s34*g44 ( j) 
g24 ( j ) =yy2-s23*g34 ( j) -s24*g44 ( j) 
gl4 ( j) =yyl-sl2*g24 ( j) -sl3*g34 ( j) -sl4*g44 ( j) 
10  continue 

do  15  j=3, jl-1 
dul=du (i, j-1, 1) 
du2=du (i, j-1, 2) 
du3=du (i, j-1, 3) 
du4=du (i, j-1, 4) 

du (i, j, 1) =du(i, j, 1) -gll ( j) *dul-gl2 ( j) 
c *du2-gl3 ( j) *du3-gl4 ( j) *du4 

du (i, j, 2) =du (i, j, 2) -g21 ( j) *dul-g22 ( j) 
c *du2-g23 ( j) *du3-g24 ( j) *du4 

du (i, j, 3) =du(i, j, 3) -g31 ( j) *dul-g32 ( j) 
c *du2-g33 ( j) *du3-g34 ( j) *du4 

du (i, j, 4) =du (i, j , 4 ) -g4 1 ( j ) *dul-g42 (j) 
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c *du2-g43 ( j) *du3-g44  ( j) *du4 

15  continue 
20  continue 
return 
end 


c 

c 

subroutine  prntff (time) 

c 

c 

common /a/ rho (100,100,2), rhou (100,100,2), rhov (100,100,2), 
c e (100,  100, 2)  , fsl (2) , fs2 (2) , 

c u (100,  100)  ,v(100,100)  ,ei(100,100)  , p (100,  100)  , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml, ggml, cv, pr, vou, ptt, ttt, pe, 

c il, jl, ilml, jlml, jlfm, n, m, kl, k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend,  nvisc 

common/b/iturb,  prt , ep (100,  100) 
common /tem/t (101, 101) , ttwall 
common / speed/ rmach (100, 100) 

write (16,10) 

10  format (lhl, 50x, 27h*  flow  field  *,//) 
write(16,20)  time, nadv 

20  format (34x, 7htime  = , el4 . 7 , 2x, 7hnadv  = ,i9) 
c 

write(50,20)  time, nadv 
c 


do  30  i=l , il, 2 
ipl=i+l 

c write (16, 24) i,  ipl 

c 24  format (lhO, 8x, 8hcolumn  , i5, 45x, 8hcolumn  ,i5) 
c write(16,25) 

c 25  format (lh  ', 3x, lhj, 5x, 3hrho, 9x, lhv, lOx,  lhu 
c c , lOx, lhp, 9x, 2hei, 12x, 

c c 3hrho, 9x, lhv, lOx, lhu, lOx, lhp, 9x,  3hei  ) 

do  29  j=l, jl 

c write  (16,  27)  j , rho  (i,  j , m)  , v (i,  j ) , u (i,  j ) , p ( i,  j ) , ei  (i,  j)  , 

c c rho (ipl, j,m) , v (ipl, j) , u (ipl, j) 

c c ,p(ipl,  j)  ,ei  (ipl,  j) 

c -insert  another  write  statement  for  graphics  file- 
c write  (50,27)  j,  rho  (i,  j,  m)  , v (i,  j)  , u (i,  j)  , p (i,  j)  , ei  (i,  j)  , 

c c rho (ipl, j,m) , v (ipl, j) , u (ipl, j) 

c c ,p(ipl,  j)  ,ei  (ipl,  j) 

c 

write  (50,  27)  j,rho(i,  j,m),v(i,j),u(i,j),p(i,j),t(i,j). 
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c rmach (i, j) , rho ( ipl , j , m)  , v ( ipl , j ) , u (ipl, j) , p (ipl,  j) , 

c t (ipl, j) , rmach (ipl, j) 

27  format(lh  , i5, 6ell . 4, 3x, 6ell . 4) 
c 

29  continue 

30  continue 
return 
end 

c 

c 

subroutine  prntxy 

c 

c 

common/ a/ rho (100,100,2), rhou (100,100,2), rhov (100,100,2)  , 
c e (100, 100,2) , fsl (2) , fs2 (2) , 

c u (100,  100) , v (100,  100) ,ei (100, 100) ,p  (100,  100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol ( 101 , 101 ) , 

c gamma, gml , ggml , cv,  pr,  vou, ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm,  n,  m, kl, k2, 

c iadbwl, eiwall, lent, cfl, 

c nadv, nend, nvisc 

ilpl=il+l 
jlpl=jl+l 

c write (16, 10) 

c 10  format  ( lhl , 50x,  38h*  mesh  coorinates  *,//) 
do  30  i=l,ilpl,4 
ipl=i+l 
ip2=i+2 
ip3=i+3 

c write(16,20)  i,  ipl,  ip2,  ip3 

c 20  format (lhO, 8x, 3 (8hcolumn  , i5, 1 6x) , 8hcolumn  ,i5,/, 
c c lh  , 4x, lhj, 6x, 3 (lhx, llx,  lhy, 16x) , lhx, llx, lhy) 

do  30  j=l, jlpl 

c write  (16,  25)  j,  z (i,  j)  , r (i,  j)  , z (ipl,  j)  , r (ipl,  j)  , z (ip2,  j)  , 

c c r (ip2,  j)  ,z  (ip  3,  j)  ,r(ip3,  j) 

25  format(lh  , i5, 3 (2el3 . 6, 3x) , 2el3 . 6) 

c -insert  another  write  statement  for  graphics  files- 

write  (50,25)  j,  z (i,  j)  , r (i,  j)  , z (ipl,  j)  , r (ipl,  j)  , z (ip2,  j)  , 
c r (ip2,  j)  , z (ip3,  j)  , r (ip3,  j) 

c 

30  continue 
return 
end 
c 

c 

subroutine  turbfl 



c 

common/a/rho (100,  100,2), rhou (100,  100,2), rhov (100,  100,2), 
c e (100,  100,2) , fsl (2) , fs2  (2) , 
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c u(100,  100) , v (100, 100) ,ei (100,  100)  , p (100,  100) , 

c fs3(2) , fs4 (2) , z (101, 101) , r (101,  101)  , vol (101,101), 

c gamma, gml, ggml, cv, pr, vou, ptt, ttt,  pe, 

c il, jl, ilml, jlml, jlfm, n,m, kl, k2, 

c iadbwl, eiwall,  lent, cf 1, 

c nadv, nend,  nvisc 

common /b/ i t urb, prt,  ep (100,  100) 
dimension  yi (100) , aomega (100) ,yyj (100) 
common /gm/thrt 

dd=2 . *thrt 
aplus=26 . 0 
ccp=l . 6 
ckleb=0 . 3 
cwk=0 . 25 
smlk=0 . 4 
bigk=0 .0168 
prt=0 . 9 
aa=0 . 80 
bb=-7 . 10e-5 

c cmutm=14.0  to  be  used  to  simulate  transition 

do  10  i=l , il 

sj=sqrt  ( (z (i  + 1,2) -z (i, 2) ) **2  + (r (i  + 1,2) -r (i,  2)  ) **2) 

szp=-  (r(i  + l,2)-r(i,2))/sj 

srp=  (z  (i+1, 2) -z (i, 2) ) /s j 

u2max=u (i, 2) *u (i, 2) +v (i, 2) *v  (i, 2) 

dudymx=2 . 0*abs (srp*u (i, 2) -szp*v (i, 2) ) *s j/vol (i, 2) 

do  5 j=2, jlml 

dy=0 . 5*abs (szp* (z (i, j+2) -z (i, j) ) +srp* (r (i,  j+2) -r (i, j) ) ) 
dudy=abs (srp* (u (i,  j + 1) -u (i, j) ) -szp* (v (i,  j + 1) -v (i, j) ) ) /dy 
if  (dudy . gt . dudymx)  dudymx=dudy 
u2=u (i, j+1) *u (i, j+1) +v(i, j+1) *v (i, j+1) 
if  (u2 . gt . u2max)  u2max=u2 
5 continue 

temp=0 . 5* (ei (i, 1) +ei (i, 2) ) /cv 
c rmu=l . 458e-6*sqrt (temp**3) / (temp+110 . 4) 

c - uf4  property 

c rmu=3 . 357e-6*sqrt (temp) / (aa  + bb*temp) 

c - UF 4 viscosity  at  2000  K 
rmu=8 . 67e-5 

ra=sqrt (rho (i, 2, 1) *rmu*dudymx) / (aplus*rmu) 
c 

yyplus=5 . 0 

if (nadv. eq. nend) then 

yy j (i) = (1 . e+3*yyplus) / (ra*aplus) 

endif 

ilm5=il-3 

if (nadv. eq. nend. and. i . eq. ilm5)  then 
do  11  j=2, jlml 

y j=abs ( sxp  * (z (ilm5, j + 1) -z (iilm5,2) ) 

* +syp* (r (ilm5, j+1) -r (ilm5, 2) ) ) 
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yplus=ra*aplus*y j 

c - calculate  emperical  uplus  using  deissler's  formula 
if  (yplus.lt. 5)  then 
uuplus=yplus 

else  if  ( yplus.gt.5  . and.  yplus.lt. 26)  then 

uuplus=-3.05  + 5 . 0*alog (yplus) 

else 

uuplus= (1./0.36) *alog (yplus ) +3 . 8 
c uuplus=5.5  + 2 . 5*alog (yplus) 

endif 

uj=0 . 5* (u (ilm5, j+1) +u (ilm5, j) ) 
uplus=u j/sqrt (rmu*dudymx/rho (ilm5, 2,1)) 
write (22,*)  yplus, uplus, uuplus 
c write(*,*)  nadv, yj , yplus, uplus 

11  continue 

endif 
fmax=0 . 0 
ymax=l . 0 
if lag=0 
do  7 j=2 , jlml 
ii=i 

if  (ii.eq.l)  ii=2 
if  (ii.eq.il)  ii=il-l 
due=0 .5* (u(ii  + l, j) -u (ii-1,  j)  ) 
dve=0 .5* (v(ii  + l, j) -v (ii-1,  j)  ) 
dze=z  (i+l,j  + l)-z(i,j  + l) 
dre=r (i+1, j+1) -r (i, j+1) 
dun=u (i, j+1) -u (i, j) 
dvn=v (i, j+1) -v (i, j) 

dzn=0.25* (z (i+1, j+2) +z (i, j+2) -z (i  + 1, j) -z  (i, j) ) 
drn=0 .25* (r (i+1, j+2) +r (i, j+2) -r (i+1, j) -r (i, j) ) 
dzr=1.0/ (dze*drn-dzn*dre) 
dvz=  (dve*drn-dvn*dre) *dzr 
dur=- (due*dzn-dun*dze) *dzr 

y j=abs (szp* (z (i, j + 1) -z (i, 2) ) +srp* (r (i, j + 1) -r  (i, 2) ) ) 

yi ( j) =yj 

aomega ( j ) =abs (dur-dvz ) 

f j=yj*aomega ( j) * (1 . 0-exp (-yj*ra)  ) 

if  (f j . gt . fmax.and. j . le. jlfm)  then 

fmax=f j 

ymax=y j 

if lag=l 

endif 

7 continue 

fwake=ymax*fmax 

if  (if lag . eq. 1 ) then 

fwakep=cwk*ymax*u2max/fmax 

if  (fwakep . It . fwake)  fwake=fwakep 

endif 

jep=jl 
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do  8 j=2 , jlml 

yj=yi  ( j) 

epi=rho (i, j, 1) * ( (smlk*yj) * (1 . 0-exp (-yj*ra)  ) ) **2*aomega ( j) 
fkleb=1.0/ (1.0+5. 5* (ckleb*y j /ymax) **6) 
epo=bigk*ccp*rho (i,  j, 1) *fwake*fkleb 
ep (i, j) =aminl (epi, epo) 
if  ( j . It . jep . and. epi . gt . epo)  jep=j 
if  (j.gt.jep)  ep  (i, j) =epo 
8 continue 
ep (i, 1) =0 . 0 
ep (i, jl) =ep(i, j 1-1 ) 

10  continue 

if (nadv . eq. nend) then 
do  100  i=2, il-2 
zd=0 .5* (z (i,2)+z (i+1, 2) ) /dd 
write(999,*)  zd,yyj(i) 

100  continue 

endif 
return 
end 
c 

c 

subroutine  rdwrt (ise, time) 

c 

c 

c subroutine  to  read  in  data  from  previous  runs  or  write 

c data, 
c 

common/a/rho (100,100,2) , rhou (100,100,2), rhov (100, 100,2) , 
c e (100, 100, 2) , fsl (2) , fs2  (2) , 

c u (100, 100) , v (100, 100) ,ei (100, 100) , p (100, 100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml , ggml , cv, pr, vou, ptt , ttt , pe, 

c il, jl, ilml, jlml, jlfm, n,m, kl, k2, 

c iadbwl, eiwall,  lent , cf 1, 

c nadv, nend, nvisc 

common/e/gll (101),gl2(101),gl3(101),gl4(101), 
c g21 (101) , g22 (101) ,g23  ( 101 ) , g24 ( 101 ) , 

c g31 (101) , g32 (101) ,g33 (101) ,g34 (101) , 

c g41 (101) , g42 (101),g43(101),g44(101), 

c ff (101) , implt 

if  (ise.eq.2)  go  to  20 

read (3, *)  nadv, time, il, jl, vou, ptt , ttt , pe 
read  (3,  *)  ( (rho  (i,j,l),u(i,j),v(i,j),ei(i,j),p(i,j), 

c vol (i, j) , i=l, il) , j=l, jl) 

ilpl=il+l 
jlpl=jl+l 

read  (3,  *)  ((z(i,j),r(i,j),  i=l,  ilpl)  , j=l,  jlpl) 

do  10  j = 1 , j 1 

do  10  i=l,il 

rho ( i , j , 2 ) =rho ( i , j , 1 ) 
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rhou  (i,  j,  1)  =rho  (i,j,l)*u(i,j) 
rhov  (i,  j,  1)  =rho  (i,j,l)*v(i,j) 
rhou ( i , j , 2 ) =rhou ( i , j , 1 ) 
rhov ( i , j , 2 ) =rhov ( i , j , 1 ) 

e (i,  j,  1)  =rho  (i,j,l)*(ei(i,j)+0.5*(u(i,j)*u(i,j) 

1 +v  (i,  j)  *v  (i,  j)  ) ) 


e(i,j,2)=e(i,j,l) 

10  continue 
return 
20  continue 

write (4 , * ) nend, time, il, jl, vou, ptt , ttt , pe 
write  (4,  *)  ( (rho  (i,j,l),u(i,j),v(i,j),ei(i,j),p(i,j), 

c vol (i, j) , i=l, il) , j=l, jl) 


ilpl=il+l 

jlpl=jl+l 

write (4,*)  ( (z (i, j) , r (i, j) , i=l, ilpl) , j=l, jlpl) 

return 

end 


c 

subroutine  fsifs (i, j, ii) 

c 

c 

common /a/ rho (100, 100,2) , rhou (100, 100,2) , rhov (100, 100,2) , 
c e (100, 100, 2) , fsl (2) , fs2 (2) , 

c u (100, 100) , v (100, 100) ,ei (100, 100) , p (100, 100) , 

c fs3 (2) , fs4 (2) ,z (101, 101) , r ( 1 0 1 , 101) , vol (101,101), 

c gamma, gml, ggml, cv, pr, vou,  ptt, ttt, pe, 

c il, jl, ilml, jlml, jlfm, n,m, kl, k2, 

c iadbwl, eiwall, lent,  cfl, 

c nadv, nend, nvisc 

common /spit /if Ixs, sz, sr 

si=sqrt (sz*sz+sr*sr) 

szp=sz/si 

srp=sr/si 

rr=rho (ii, j,  n) 

uu=u (ii, j) 

vv=v  (ii, j) 

up=  szp*uu+srp*vv 

vp=-srp*uu+szp*vv 

cc=sqrt (ggml*abs (ei  (ii, j) ) ) 

alp=0 . 5* (up*up+vp* vp) 

do  10  nn=l,2 

il=i+2-nn 

i2=il+3-2*nn 

if  (if lxs . eq. 1 . or . i . le . nn . or . i . ge . il+nn-3)  i2=il 
if  (abs  (p(i2,  j)-p(il,  j)  ) .gt. 
c 0.5*aminl (abs (p (il, j) ) , abs (p (il, j) ) ) ) i2=il 
if  (abs (p (i  + 1, j) -p  (i, j) ) .gt . 
c 0.5*aminl (abs (p (i+1, j) ) , abs (p (i, j) ) ) ) then 
i2=il 
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rr=rho (il, j, n) 

uu=u (il, j) 

vv=v (il, j) 

up=  szp*uu+srp*vv 

vp=-srp*uu+szp*vv 

cc=sqrt (ggml*abs (ei (il,  j) ) ) 

alp=0 . 5* (up*up+vp*vp) 

endif 

dr=l . 5*rho (il, j,  n) -0 . 5*rho (i2, j, n) 
dru=l . 5*rhou (il, j, n) -0 . 5*rhou (i2, j, n) 
drv=l . 5*rhov (il, j, n) -0 . 5*rhov (i2, j, n) 
de=l . 5*e (il, j,n)-0.5*e(i2,  j,  n) 
dp= (dr*alp-uu*dru-vv*drv+de) *gml 

xxl=dr-dp/ (cc*cc) 

xx2= (-up*dr+szp*dru+srp*drv) *cc+dp 
xx3= (-vp*dr-srp*dru+szp*drv) /rr 
xx4=-xx2+2 . 0*dp 


ddl=si*up 
dd2=si* (up+cc) 
dd4=si* (up-cc) 

sgns=2*nn-3 

yyl=0 . 5* (ddl+sgns*abs (ddl) ) *xxl 
yy2=0 .5* (dd2+sgns*abs (dd2) ) *xx2 
yy3=0 . 5* (ddl+sgns*abs (ddl ) ) *xx3 
yy4=0 . 5* (dd4+sgns*abs (dd4 ) ) *xx4 

t0=yyl+0.5* (yy2+yy4)  / (cc*cc) 

t l=up*t 0+0 . 5* (yy2-yy4) /cc 

t2=vp*t0+rr*yy3 

fsl  (k2 ) = (nn-1 ) *fsl (k2)  +t0 

fs2 (k2) = (nn-1) *fs2 (k2) +szp*tl-srp*t2 

fs3  (k2)  = (nn-1) *fs3(k2)+srp*tl  + szp*t2 

fs4 (k2) = (nn-1) *fs4 (k2) + 

c alp*t0+rr*vp*yy3+0 .5* (up* (yy2-yy4) /cc+ (yy2+yy4) /gml) 

10  continue 
return 
end 


subroutine  fs jfs (i, j, j j) 

c 

c 

common/a/rho (100, 100,2), rhou (100, 100,2), rhov (100, 100,2), 
c e (100,  100,2) , fsl (2) , fs2  (2) , 

c u (100, 100) , v (100, 100) ,ei(100,100) , p (100, 100) , 

c fs3 (2) , fs4 (2) , z (101, 101) , r (101, 101) , vol (101,101), 

c gamma, gml, ggml, cv,pr,vou,ptt,ttt,pe, 

c il, jl, ilml, jlml, jlfm, n,m, kl, k2, 

c iadbwl, eiwall, lent, cfl, 


c nadv, nend, nvisc 

common/splt /if lxs, sz, sr 
if  ( j .eq. 1 .or. j .eq. jlml)  return 
sj=sqrt (sz*sz+sr*sr) 
szp=sz/s j 
srp=sr/s j 

rr=rho (i, j j,n) 
uu=u  (i,  j j) 
vv=v ( i , j j) 

cc=sqrt (ggml*abs (ei (i, j j) ) ) 

do  10  nn=l,2 

up=  srp*uu-szp*vv 
vp=  szp*uu+srp*vv 
alp=0 . 5* (up*up+vp*vp) 


jl=j+2-nn 
j2= jl+3-2*nn 

if  (if lxs . eq. 1 . or . j . le . nn. or . j . ge . jl+nn-3)  j 2 = j 1 
if  (abs  (p(i,  j2)-p(i,  jl)  ) .gt. 
c 0 . 5*aminl (abs (p (i, jl) ) , abs (p (i, j2) ) ) ) j 2= j 1 
if  (abs (p (i, j + 1) -p  (i, j) ) .gt . 
c 0.5*aminl (abs (p (i, j+1) ) , abs (p (i, j) ) ) ) then 
j 2= j 1 

rr=rho (i, jl,  n) 

uu=u (i, jl) 

vv=v (i, jl) 

up=  srp*uu-szp*vv 

vp=  szp*uu+srp*vv 

cc=sqrt (ggml*abs (ei (i, jl) ) ) 

alp=0 . 5* (up*up+vp*vp) 

endif 

dr=l . 5*rho (i, jl, n) -0 . 5*rho (i, j2,  n) 
dru=l . 5*rhou  (i,  jl,  n)  -0 . 5*rhou  (i,  j2,  n) 
drv=l . 5*rhov (i, jl , n) -0 . 5*rhov (i,  j2 , n) 
de=l . 5*e  (i, jl,n)-0.5*e(i,j2,n) 
dp= (alp*dr-uu*dru-vv*drv+de) *gml 

wtl=l . 0 
wt2=0 . 0 

if  (j.lt.jlfm)  then 

wtl=(  (z  (i,  j)-z  (i,2)  ) **2+(r  (i,  j)-r  (i,2)  ) **2)  / 
c ( (z (i,  jlfm) -z (i,2) ) **2+ (r (i,  jlfm) -r (i, 2) ) **2) 
wt2=l . 0-wt 1 

dr=wtl*dr+wt2*rho (i, jl, n) 
dru=wtl*dru+wt2*rhou (i, jl,  n) 
drv=wtl*drv+wt2*rhov (i, jl,n) 
de=wtl*de+wt2*e (i, jl,  n) 
dp=wtl*dp+wt2*p (i,  jl) 
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endif 

xxl=dr-dp/ (cc*cc) 
xx2= (-up*dr+srp*dru-szp*drv) /rr 
xx3= (-vp*dr+szp*dru+syr*drv) *cc+dp 
xx4=-xx3+2 . 0*dp 

ddl=s j*vp 
dd3=s j* (vp+cc) 
dd4=s j* (vp-cc) 

sgns=2*nn-3 

yyl=0.5* (ddl+sgns*abs (ddl) ) *xxl 
yy2=0 . 5* (ddl+sgns*abs (ddl) ) *xx2 
yy3=0 .5* (dd3+sgns*abs (dd3) ) *xx3 
yy4=0 . 5* (dd4+sgns*abs (dd4 ) ) *xx4 

vpj=szp*u(i, jl)+srp*v(i,  jl) 
vp=wt l*vp+wt2*vp j 
alp=0.5* (up*up+vp*vp) 

t0=yyl+0.5* (yy3+yy4) / (cc*cc) 

tl=up*tO+rr*yy2 

t2=vp*t0+0 .5* (yy3-yy4 ) /cc 

fsl (k2) = (nn-1) *fsl (k2) +tO 

fs2 (k2) = (nn-1) *fs2 (k2) +srp*t l+szp*t2 

fs3  (k2)  = (nn-1) *fs3 (k2) -szp*tl  + srp*t2 

fs4 (k2)  = (nn-1 ) *fs4  (k2)  + 

c alp*t0+rr*up*yy2+0 . 5* (vp* (yy3-yy4) /cc+ (yy3+yy4) /gml) 
10  continue 
return 
end 
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